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ABSTRACT 

A  numerical  technique  is  formulated,  in  a  computer  program  U2DIIF,  for  the 
solution  of  flow  over  an  airfoil  executing  an  arbitrary  unsteady  motion  in  an  inviscid 
and  incompressible  medium.  The  technique  extends  the  well  known  Panel  Methods  for 
steady  flow  into  solving  a  non-linear  unsteady  flow  problem  arising  from  the 
continuous  vortex  shedding  into  the  trailing  wake  due  to  the  unsteady  motion  of  the 
airfoil.  Numerous  case-runs  are  presented  to  verify  U2DIIF  computer  code  against 
other  theoretical  and/ or  numerical  methods  as  well  as  in  cases  where  limited 
experimental  data  are  obtainable  in  literatures.  These  case-runs  include  airfoils 
undergoing  a  step  change  or  a  modified  ramp  change  of  angle-of-attack,  airfoils 
executing  harmonic  oscillation  in  pitching  and  plunging  motions  and  airfoils 
penetrating  a  sharp  edge  gust. 
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THESIS  DISCLAIMER 

The  reader  is  cautioned  that  computer  programs  developed  in  this  research  may 
not  have  been  exercised  for  all  cases  of  interest.  While  every  effort  has  been  made, 
within  the  time  available,  to  ensure  that  the  programs  are  free  of  computational  and 
logic  errors,  they  cannot  be  considered  validated.  Any  application  of  these  programs 
without  additional  verification  is  at  the  risk  of  the  user. 
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I.  INTRODUCTION 

A.  GENERAL 

In  this  thesis,  a  numerical  method  is  formulated  and  coded  in  a  FORTRAN 
computer  program,  codename  U2DIIF  (Unsteady  2-Dimensional  Inviscid 
Incompressible  Flow),  to  solve  for  the  flow  over  an  airfoil  which  is  executing  an 
unsteady  time-dependent  motion  in  an  inviscid,  incompressible  medium. 

B.  APPROACH 

The  basic  approach  to  this  problem  is  the  extension  of  a  very  general  and 
powerful  technique,  called  Panel  Methods,  developed  by  Hess  &  Smith  [Ref.  1]  for 
steady  potential  flow  problems,  to  include  the  unsteady  motion  of  the  airfoil  that  is 
continuously  shedding  vorticity  into  the  trailing  wake.  This  vortex  shedding  process 
creates  the  non-linearity  effects  of  the  problem  in  that  the  wake  vortices  influence  the 
flow  over  the  airfoil  which  in  turn  alters  the  vortex  shedding  as  the  airfoil  proceeds  in 
time.  It  is  this  very  non-linearity  of  unsteady  flow  that  distinguishes  itself  from  the  well 
known  steady  Panel  Methods  solution  where  the  mathematical  formulation  of  the 
problem  results  in  a  set  of  N  linear  equations  in  N  unknowns  which  are  solved  easily 
with  the  standard  Gaussian  elimination  algorithm. 

The  unsteady  flow  problem  is,  however,  deprived  of  this  relatively  easy  solution 
technique.  Instead,  an  iterative  type  of  solution  is  needed  for  this  non-linear  problem. 
The  correct  mathematical  model  must  therefore  be  formulated  to  describe  the  vortex 
shedding  process  that  provides  the  mechanism  for  the  iteration  to  proceed  towards  a 
converged  set  of  solution  in  each  time  step. 

It  is  the  objective  of  this  thesis  to  develop  a  numerical  computer  program  that 
performs  this  non-linear  potential  flow  calculation  which  proceeds  step  by  step  in  time. 
At  each  time  step,  a  complete  set  of  potential  flow  solutions,  inclusive  of  the    lirfoil 
pressure   distribution,   force   and   moment   coefficients,   and   the   trailing   vortex    \ 
pattern  (strengths  and  positions  of  shed  vortices),  is  obtained. 

C.  SCOPE 

The  Panel  Method  of  Hess  &  Smith,  which  utilises  both  the  distributed  sources 
and  vorticities  as  panel  singularities,  for  steady  flow  solution  is  described  in  Chapter  II. 
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Chapter  III  formulates  the  mathematical  model  for  the  unsteady  flow  problem 
and  its  solution  procedures,  highlighting  the  essential  features  in  solving  the  non-linear 
problem  of  unsteady  flow. 

Chapter  IV  describes  the  computer  program  U2DIIF,  its  essential  capabilities, 
limitations  and  the  necessary  input  set-up  for  typical  case-runs. 

The  results  of  some  of  the  case-runs  are  presented  in  Chapter  V.  They  are 
compared  with  other  theoretical  and/ or  numerical  methods  as  well  as  in  cases  where 
limited  experimental  data  are  obtainable  in  the  literature.  These  case-runs  include 
airfoils  undergoing  a  step  change  or  a  modified  ramp  change  in  angle-of-attack,  airfoils 
executing  harmonic  oscillation  in  both  pitching  and  plunging  motions  and  airfoils 
penetrating  a  sharp  edge  gust. 

In  the  concluding  remarks  of  Chapter  VI,  the  future  development  and  application 
potential  of  this  numerical  method  to  other  studies  of  unsteady  2-dimensional  inviscid 
incompressible  flow  are  mentioned. 
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H.  STEADY  FLOW  PROBLEM  FORMULATION 

A.  FRAME  OF  REFERENCE 

Consider  a  2-dimensional  airfoil  in  motion  with  constant  linear  velocity  —  Vqq  as 
shown  in  Figure  2.1.  Using  an  (x,y)  coordinate  system  fixed  on  the  airfoil,  where  the  x- 
axis  coincides  with  the  chord  line  originating  from  the  leading  edge  towards  the  trailing 
edge  of  the  airfoil,  the  flow  in  this  frame  of  reference  is  steady.  That  is  to  say,  the  fluid 
velocity  and  pressure  in  the  flow  field  depend  only  on  the  spatial  coordinates  (x,y)  and 
not  on  time.  The  airfoil  then  appears  to  be  submerged  in  an  onset  flow  whose  velocity 
is  Vqq  and  making  an  angle  of  attack,  a,  with  the  x-axis  (see  Figure  2.1). 

B.  STEADY  FLOW  PANEL  METHODS 

1 .  Definition  of  Nodes  and  Panels 

The  airfoil  surface  is  divided  into  (n)  straight-line  segments,  called  panels,  by 
(n+  1)  arbitrary  chosen  points,  called  nodes,  distributed  over  the  airfoil  contour  as 
shown  in  Figure  2.2.  The  panel  numbering  sequence  starts  with  panel  1  on  the  lower 
surface  at  the  airfoil  trailing  edge  and  proceeds  clockwise  around  the  airfoil  contour  so 
that  the  last  panel  (panel  n)  ends  on  the  upper  surface,  also  at  the  airfoil  trailing  edge. 

Notice  that  this  numbering  sequence  dictates  that  the  airfoil  body  always  lies 
on  the  right  hand  side  of  the  1th  panel  as  one  proceeds  from  the  Ith  node  to  the  (i+  0th 
node.  Also  the  1st  and  the  (n+1)111  nodes  coincide  at  the  trailing  edge.  It  therefore 
facilitates,  as  shown  in  Figure  2.2,  the  common  definition  of  unit  normal  vector  n.  and 
the  unit  tangential  vector  t.  for  all  panels,  i.e.,  n.  is  directed  outward  from  the  body  into 
the  flow  and  t  is  directed  from  the  Ith  node  to  the  (i+  0th  node. 

2.  Distribution  of  Singularities 

Figure  2.2  also  indicates  that  a  uniform  source  distribution  q.  and  a  uniform 
vorticity  distribution  y  are  placed  on  the  jth  panel.  The  source  strength  q.  varies  from 
panel  ;o  panei  -.vtiereas  the  voracity  strength  y  remains  the  same  for  ail  panels.  I 
particular  choice  of  singularity  distributions  is  one  of  the  many  types  of  singularity 
combinations  (it  happened  to  be  the  pioneering  one  though)  ever  used  in  a  wide  variety 
of  the  so  called  Panel  Methods.  The  success  of  representing  the  flow  past  an  arbitrary 
shaped  airfoil  by  surface  singularity  distributions  lies  in  the  fact  that  these  singularity 
distributions  automatically  satisfy  Laplace's  equation,  the  governing  flow  equation  for 
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Figure  2. 1     Frame  of  Reference  for  Steady  Flow. 
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y 


i+1 


n.  =   —  sinG.  i  +  cosGj  j 
t.  =     cosG.  i  +  sinG.  j 


i+i  t 


Panel  j 
-<s  Source  Distribution  q. 
Vorticity  Distribution  y 


Figure  2.2     Panel  Methods  Representation  for  Steady  Flow. 
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inviscid  incompressible  flow,  and  the  boundary  condition  at  the  far  field  (oo).  In 
addition,  the  superposition  principle  applies  to  any  linear  homogeneous  second  order 
partial  differential  equation  such  as  Lapiace's  equation.  Therefore  one  can  build  up  an 
overall  complicated  flow  field  by  the  combination  of  simple  flows  if  the  appropriate 
boundarv  conditions  on  the  airfoil  can  be  satisfied  accuratelv.  For  our  case  the  overall 
flow  field  (represented  by  the  velocity  potential  <P)  can  be  built  up  by  three  simple 
flows, 

<D  =  <Poo  +  <PS  +  <PV  (eqn2.1) 

where  (Poo  is  the  potential  of  the  onset  flow, 

(Pqq  =  Vqq  (x  cosa  +  y  sina)  (eqn  2.2) 

<ps  is  the  velocity  potential  of  a  source  distribution  of  strength  q(s)  per  unit  length, 

<ps  =  J  — In  r  ds  (eqn  2.3) 

cp   is  the  velocity  potential  of  a  vorticity  distribution  of  strength  y(s)  per  unit  length. 


<ps  =   -  J  ^——  6  ds  (eqn  2.4) 


The  integrals  in  Equations  2.3  and  2.4  are  performed  along  the  surface 
contour  s  and  (r,0)  are  polar  coordinates  of  any  field  point  (x,y)  measured  from  the 
airfoil  surface  at  an  arbitrary  point  as  shown  in  Figure  2.3.  The  difficult  task,  of 
evaluating  these  integrals  has  been  greatly  simplified  by  our  singularity  distributions 

postulated  to  represent  the  flow  over  the  airfoil:  that  is,  instead  of  integrating  over  the 
entire  airfoil  contour,  we  integrate  on  each  panel  along  a  straight  line  where  q.  and  y 
are  constant,  then  sum  up  the  effects  of  all  panels.  Equation  2.1  therefore  becomes, 

O  =  Vqo  (x  cosa  +  y  sina)  +  £  Jpane, .  [  JjJL  in  r  -  -~-  6  ]  ds       (eqn  2.5) 
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Field  Point 

(x,y) 


O  =  V  —  (x  cosa  +  y  sina)    + 


J  — In  r  ds   -   f G  ds 

J    2k  j   2k 


Figure  2.3     Potential  Evaluation  at  a  Field  Point. 


20 


It  can  be  seen  from  Equation  2.5  that  <l>  is  completely  defined  if  the  (n+  1) 
unknowns,  q.  (j=  1,2,. ...,n)  and  y,  can  be  calculated  using  a  numerical  technique  yet  to 
be  described.  Once  the  potential  O  is  solved,  the  velocity  can  be  evaluated  by  taking 
grad<&.  At  this  point  we  introduce  a  definition  of  disturbance  potential,  <p,  as  the  sum 

of  potential  due  to  both  the  source  and  vorticity  distribution, 

<p  =  (ps  +  (pv  (eqn  2.6) 

Equation  2.1  therefore  reads, 

<&  =  <Poo  +  <P  (e3n  2J) 

The  total  velocity  vector  is  thus, 

total 

=  ^(Pco  +  Vq> 

=  voo  +  V(P  (ecin  2-8) 

The  pressure  can  be  obtained  from  Bernoulli's  Equation, 


P-P™  V      ,    O 

C    =  °°  =  1  -  ( -^  )2  (eqn  2.9) 


'P        'ApVoo2  V 


oo 


Notice  that  Figure  2.3  indicates  that  the  field  point  lies  off  the  airfoil  surface, 

however,  we  are  interested  in  field  points  that  are  on  the  airfoil  surface.  In  the  case  of 

steady  flow,  the  expressions  for  VtoLal  and  C    are  the  same  for  field  points  lying  on  or 

off  the  airfoil  surface.  It  is  nevertheless  not  the  same  in  unsteady  flow,  as  will  be  seen 

in  Chanter  III,  in  that  Vf     ,  must  include  the  rigid  bodv  motion  of  the  airfoil  when  one 

total 

evaluates  field  points  on  the  airfoil  surface. 
3.  Boundary  Conditions 

The  boundary  conditions  to  be  satisfied  include  the  flow  tangency  conditions 
and  the  Kutta  Condition.  The  flow  tangency  conditions  are  satisfied  at  the  exterior  mid 
points,  called  control  points,  of  all  panels  by  taking  the  resultant  velocity  at  each 
control  point  to  have  only  (V1).  but, 
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(V11).  =  0  ,    i=l,2 ,n  (eqn2.10) 

where  (V1).  and  (Vn).  are  the  tangential  and  normal  components  of  the  total  velocity  at 
the  control  point  of  the  Ith  panel  due  to  the  free  stream  and  the  velocities  induced  by 

the  source  and  vorticity  distributions  on  all  the  panels,  j  (j=  1,2, n). 

The  Kutta  condition  postulates  that  the  pressures  on  the  upper  and  lower 
panels  at  the  trailing  edge  be  equal  in  order  that  the  flow  leaves  the  trailing  edge 
smoothly.  By  using  Bernoulli's  equation  for  steady  potential  flow,  this  pressure 
equilibrium  condition  implies  that  the  tangential  velocities  in  the  downstream  direction 
at  the  1st  and  the  n*  panel  control  points  must  be  equal.  This  fact  is  certainly 
consistent  with  the  knowledge  that  when  steady  flow  is  established,  the  total  circulation 
over  the  airfoil  does  not  change  if  the  tangential  velocities  are  the  same  at  the  trailing 
edge  panels. 

(V^  =    -(Vl)n  (eqn2.11) 

If  one  could  explicitly  express  Equations  2.10  and  2.11  in  terms  of  the  unknowns 
q.  (j=  l,2,....,n)  and  y,  the  task  is  then  reduced  to  solving  a  linear  system  of  (n-M) 
simultaneous  equations  for  the  (n+  1)  unknowns. 

C.       INFLUENCE  COEFFICIENTS 

1.  The  Concept  of  Influence  Coefficients 

The  numerical  technique  employed  in  Panel  Methods  to  manipulate 
equations  2.10  and  2.11  into  an  algebraic  system  of  linear  simultaneous  equations 
involves  the  important  concept  of  influence  coefficients.  An  influence  coefficient  is 
defined  as  the  velocity  induced  at  a  field  point  by  a  unit  strength  singularity  (be  it  a 
point  singularity  or  a  distributed  singularity)  placed  anywhere  within  the  flow  field.  In 
this  case,  it  is  the  unit  strength  singularity  distribution  on  one  panel.  Recall  that 
equations  2.10  and  2.11  simply  require  the  computation  of  the  normal  and  tangential 
velocity  components  at  all  the  panel  control  points.  The  normal  components  oi' 
velocities  are  essential  in  satisfying  How  tangency  conditions  while  the  tangential 
components  of  velocities  are  necessary  for  satisfying  the  Kutta  condition  as  well  as 
computing  the  pressure  distribution.  The  procedure  is  thus  to  compute,  at  the  i^  panel 
control  point,  the  velocity  components  induced  by  the  source  and  vorticity 
distributions  on  all  the  panels,  j  (j  =  1,2 n),  including  the  Ith  panel  itself.  Summation 


22 


of  all  the  induced  velocities,  separately  for  the  normal  and  tangential  components, 
together  with  the  free  stream  velocity  components  produces  all  the  required  (Vn).  and 
(V%  i-  1,2 ,n. 

2.  Notation  for  Influence  Coefficient 

We  shall  adopt  a  consistent  set  of  notation  for  the  influence  coefficients  used 
throughout  this  documentation.  It  is  so  designated  to  permit  easy  recognition  in  that 
each  influence  coefficient  contains  all  the  associated  information  one  needs.  An 
influence  coefficient  is  denoted  with  a  superscript  and  two  subsript  as  follows: 

k  pq 

where  x  denotes  the  type  of  singularity  involved,  we  shall  arbitrarily  use  A,  B  and  C  for 

the  uniformly  distributed  source,  uniformly  distributed  vorticity  and  point  vortex 
respectively.  The  superscript  s  is  an  indicator  ceiling  which  component  the  induced 
velocity  is.  The  first  subscript  p  identifies  the  field  point  where  the  induced  velocity  is 
evaluated.  The  second  subscript  q  denotes  the  particular  singularity  contributing  to  the 
induced  velocity. 

We  thus  define,  for  the  steady  flow  problem,  the  following  influence 
coefficients  : 

•  Anj.  :  normal  velocity  component  induced  at  the  i*  panel  control  point  by  unit 
strength  source  distribution  on  the  Ith  panel. 

•  Alr  :  tangential  velocity  component  induced  at  the  Ith  panel  control  point  by 
unit  strength  source  distribution  on  the  j™  panel. 

•  B".  :  normal  velocity  component  induced  at  the  i*  panel  control  point  by  unit 
strength  vorticity  distribution  on  the  j™  panel. 

•  B1..  :  tangential  velocity  component  induced  at  the  Ith  panel  control  point  by 
unit  strength  vorticity  distribution  on  the  j1^  panel. 

3.  Computation  of  Influence  Coefficients 

The  influence  coefficients  turn  out  to  be  related,  not  surprisingly,  to  the 
geometry  of  the  airfoil  and  the  manner  in  which  the  panels  are  formed.  Specifically,  as 
derived  in  [Ref.  2],  the  A's  and  B's  influence  coefficients/  due  to  uniformly  distributed 
source  or  vorticity  are  functions  of : 

•  The  natural  logarithm  of  the  ratio  of  distance  from  the  i*  panel  control  point 
(the  field  point)  to  the  (j+  l)^  and  j"1  nodes  of  the  Ith  panel  where  singularities 
are  distributed. 


C's  coefficients  will  be  needed  only  for  unsteady  flow. 
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•  The  angle,  in  radian,  subtended  at  the  Ith  panel  control  point  (the  field  point)  by 
the  (j  -  l)th  and  j     nodes  of  the  j     panel  where  singularities  are  distributed. 

•  The  trigonometry  angles  of  the  Ith  and  j"1  panels. 

Referring    to    the    geometrical    quantities    indicated    in     Figure    2.4,    the 
expressions2  for  these  influence  coefficients  are  : 


27T 

1 


[  sin(9.  -  G.)  In  -il±i.  +  cos(9.  -  6.)  p..  ]  ,        for  i  *  j 

,        for  i  =  j      (eqn  2.12) 


^t      —   .      __  r  *i~(Q   _  Q  \  R       —    „~cft\   -Q\U  __AJ_±JL 


..  =  [  sin(6.  -  6.)  p..  -  cos(0.  -  0.)  In     U'^L  ]  ,        for  i  x  j 

.j        2-k  x       ]     l]  l       ]  r.. 

=  0  ,        for  i  =  j       (eqn  2.13) 


1  r.  .  .  , 

Bn.  =  [cos(9.-0.)  In     '''~!    -  sini'9. -8.)  p..  '  .        for  i  x  j 

i]        2k    '  l       J  r..  l       J      :J 

=  0  for  i  =  j       (eqn  2.14) 

« 

B1..  -  [  cos(0.  -  0.)  p..  +  sin(0.  -  0.)  In  ili+L  j  for  i  x  j 

i]         27i  l       J     *1  l       V  r.. 

,        for  i  =  j        (eqn  2.15) 


1 

T 


where  : 


ry+i  =  V[(xml-xj  +  1)2  +  (yxnl-yj  +  1)2] 

r.j  =  V  [(xmj-xp2  +  (ymj-Vj)2] 

ym,  =  K(yi+yl+1) 


Actual  computation  uses  A1-.  =    -  Bn..   and   B1..  =  A"-,   to  reduce  computing 

time. 
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Figure  2.4     Influence  Coefficients  due  to  Uniformly  Distributed  Singularities. 
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9.  =  arciant— — ! — -L~) 


x.  4_ ,  —  x, 


8.  =  arctan(  y'  +  1     ^ 


x.  _i_ ,  —  x. 

j  +  i       j 


■) 


ym.  -  y.  .  ,  ym.  —  y. 

p..  =  arctan(       1       '  +  1)  -  arctan(       '       ') 


xmi-xj  +  1  xm[-x] 


D.       NUMERICAL  SOLUTION  SCHEME 

1.  Rewriting  the  Boundary  Conditions 

Using  the  concept  of  influence  coefficients,  the  flow  tangency  conditions  of 
Equation  2.10  can  be  expressed  as, 

n  n 

I  [  An;j  q,  ]  -  Y  I  Bnj.  +  V^o  sin(a  -  9.)  =  0  ,        i-  l,2,....,n  (eqn  2.16) 

i  - 1  i  - 1 

The  Kutta  condition  of  Equation  2.11,  in  terms  of  influence  coefficients,  looks 
as, 

n  n 

j=i  j=i 

n  n 

=  I  [  Acnj  qj  ]  +  Y  Z  Bcnj  +  Voo  cos(a-6n)  (eqn  2.17) 

l-l  i-l 

The  negative  signs  appearing  on  the  left-hand-side  of  Equation  2.17  are  a 
direct  consequence  of  our  definition  of  unit  tangential  vector.  In  other  words,  the 
tangential  velocities  on  the  lower  surface  panels  downstream  of  the  front  stagnation 
point  have  negative  values,  rhis  feature  in  fact  allows  one  to  predict  the  front 
stagnation  point  by  interpolating  the  velocity  distribution  around  the  leading  edge. 

2.  Solving  the  Strengths  of  Source  and  Vorticity  Distributions 

It  is  not  difficult  at  this  stage  to  see  that  if  we  collect  the  like  terms  in 

Equation   2.17  and  expand   Equation   2.16  for  all  i's  (i=l,2, n),   these   equations 

constitute  none  other  than  a  linear  algebraic  system  of  (n  +  1)  equations  as  shown  in 
the  matrix  Equation  2.18. 
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al,l       al,2       al,3       aU  +  l 

a2,l       a2,2       a2,3       a2,n  +  l 

a3,l       a3,2       a3,3       ••••;•    a3,n+l 

an,l       an,2        an,n+l 

an  +  l,l       an  +  l,n+l 


*1 

^2 
^3 


b, 

b, 


;n+l 


(eqn2.18) 


Equation  2.18  is  a  set  of  linearly  independent  equations  which  can  be  easily 
solved  by  any  standard  linear  system  solver,  one  of  which  is  the  well  known  method  of 
Gaussian  Elimination  with  Partial  Pivoting. 

3.  Computation  of  Velocity  and  Pressure  Distribution 

Once  the  q.  (j=  l,2,....,n)  and  y  are  solved,  the  velocities  at  all  the  panel 
control  points  can  be  evaluated.  Only  the  tangential  components  exist  since  the  normal 
components  are  already  set  to  zeroes  by  the  flow  tangency  conditions. 


total 


V 


CO 


—     f\Tt 


-  (Vl)j   ,        i=l,2,....,n 


(eqn2.19) 


where 


n 


n 


(Vl)j  =  £  [  A\.  q.  ]  +  y  I  B\.  +  Voo  cos(a  -  8.)   ,        i=  l,2,....n         (eqn  2.20) 

Substituting  Equation  2.20  into  the  C    equation  (Equation  2.9),  the  pressure 
coefficients  at  the  Vth  panel  control  point  is  : 


(Cp).  =  1  -  (V1).2   ,        i-  1.2. n 


(eqn  2.21) 


4.  Computation  of  Forces  and  Moments 

The  two   dimensional   aerodynamic   coefficients   of  lift  (Cp),   drag  (Cd)   and 
pitching  moment  (C   )  about  the  leading  edge  are  computed  by  integration  of  the 


pressure  distribution  assuming  constant  C    exists  in  each  panel.  The  computation  is 
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first  done  by  integrating  forces  in  the  airfoil-fixed  coordinate  system  followed  by  a 
rotation  to  the  respective  lift  and  drag  directions  along  and  perpendicular  to  the  free 
stream  (V-jq)  as  follows  : 


S"£(C|Afr|  +  l-*i)  (eqn2.22) 

i=l 
Cx  =    -I^i^  +  l"^  (eqn2.23) 

i—  1 
n 

cm  =  I  <cP)i  «xi + 1  -  *i)  xmi  +  fri + 1  -  y-)  ymi  ]  (*w 2-24) 

i=l 


Cd  =  Cx  cosa  +  C    sina  (eqn  2.25) 


Cp  =  C„  cosa  -  Cv  sina  (eqn  2.26) 

v  y  x 
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III.  UNSTEADY  FLOW  PROBLEM  FORMULATION 

A.       OVERVIEW  OF  UNSTEADY  FLOW  MODELING 

1.  Some  Previews 

Having  fully  understood  the  Panel  Methods  formulation  and  solution  for  the 
steady  flow  problem,  one  could  then  venture  into  the  interesting  and  complicated 
unsteady  flow  case.  In  this  Chapter,  we  shall  see  how  we  could  build  the  time- 
dependency  into  the  Panel  Methods  solution  which  has  been  proven  to  be  an  useful 
and  accurate  tool  for  steady  flow.  The  approach  in  the  unsteady  flow  problem 
formulation  will  proceed,  in  general,  in  a  manner  similar  to  Chapter  II.  However,  as  we 
go  along,  we  will  pick  up  the  highlights  of  the  essential  differences  (also  similarity) 
between  the  two  problems.  Additional  flow  modeling  of  the  vortex  shedding  process 
that  greatly  influences  the  numerical  solution  technique    will  be  discussed  in  details. 

2.  Specific  Unsteady  Flow  Model 

Recall  that  in  steady  flow,  the  problem  is  considered  solved  as  soon  as  the 
airfoil  surface  singularity  distributions  of  source  and  vorticity  q.  (j=  l,2,....,n)  and  y  are 
determined.  These  (n  +  1 )  unknowns  are,  however,  time  dependent  in  unsteady  flow. 
We  therefore  introduce  a  subscript  k  as  the  time-step  counter;  that  is,  we  postulate  to 
solve  the  unsteady  flow  problem  at  successive  intervals  of  time,  starting  from  tQ  =  0.  At 
each  time-step  t.  (k=  1,2,...., 00),  we  represent  the  airfoil  by  surface  singularity 
distributions  consisting  of  source  distribution  (q.).  (j=  l,2,....,n)  and  vorticity 
distribution  yk>  Again  the  source  strengths  vary  from  panel  to  panel  but  the  vorticity 
strength  remains  the  same  for  all  panels. 

The  overall  circulation  Tk  at  time-step  tfe  is  simply  y.  multiplied  by  the  airfoil 
perimeter,  I.  Since  the  total  circulation  in  a  potential  flow  field  must  be  preserved 
according  to  the  Helmholtz's  theorem  of  continuity  of  vorticity,  any  changes  in  the 
circulation  on  the  airfoil  surface  must  be  manifested  by  an  equal  and  opposite  change 
in  vorticity  in  "he  wake.  We  call  this  the  vortex  shedding  process  and  postulate,  as 
shown  in  Figure  3.1,  that  this  shed  vorticity  takes  place  through  a  small  straight  line 
wake  element  attached  as  an  additional  panel  to  the  trailing  edge  with  uniform  vorticity 
distribution  (yw)k>  We  shall  from  now  on  refer  to  this  panel  as  the  shed  vorticity  panel, 
The  shed  vorticity  panel  will  be  established  if  its  length  Ak  and  inclination  0k,  to  x- 


Referring  to  the  switch  from  a  direct  scheme  to  an  iterative  scheme. 
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Vortex  Shedding  at  Time  Step  tk 


Panel  j 

<   Source  Distribution  (q.)k 
Vorticity  Distribution  y. 


j.  i')+[ 


Heimholtz's  theorem 


\  (yA  +  rk  =  r 
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tanO,  =  i-sk 
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Figure  3.1     Extension  of  Panel  Methods  Representation  for  Unsteady  Flow. 
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axis  of  the  airfoil-fixed  coordinate  system,  satisfy  the  Helmhoitz's  theorem  as  follows, 

\  CTA  +  Tk  =  rk.,  (eqn3.1) 

°r        \  (VA  =  rH  -  Tk  =  e  (yw  -  Yk)  (eqn  3.2) 

« 

where  Tk_1  and  y.  ,  are  respectively  the  total  circulation  and  vorticity  strengths  which 
are  already  determined  at  a  time-step  L  ,  before  t . 

Let  us  project  one  time  step  ahead  to  t.+  ,,  we  allow  the  shed  vorticity  panel 
to  be  detached  from  the  trailing  edge  and  get  convected  downstream  as  a  concentrated 
free  vortex,  with  circulation  Ak  (Yw)k  or  Fk_1  -  r"k,  according  to  the  resultant  local 
velocity  occurred  at  the  center  of  the  vortex  core.  At  the  same  time  a  brand  new  shed 
vorticity  panel  is  formed  for  the  new  time  step  and  the  process  is  repeated.  Therefore 
the  shed  vorticity  panel  model  provides  exactly  the  desired  communication  mechanism 
to  carry  the  solution  from  one  time-step  to  another. 

We  now  return  back  to  the  time- step  t.  and  immediately  realise  that  as  a 
result  of  this  continuous  vortex  shedding,  there  has  been  a  series  of  shedding  processes 
occurred  prior  to  tk  that  cummulated  in  a  string  of  concentrated  core  vortices  of 

strengths  (r"k_2  -  rfc-1),  (^.3  _  rk.2)>  (rk.4-rk.3),    and  so  on,  forming  the 

wake  pattern  behind  the  airfoil  as  shown  in  Figure  3.1 

The  presence  of  the  shed  vorticity  panel  and  the  downstream  resultant  wake 

core   vortices   do   influence   the   upstream  flow  in   inviscid   incompressible   flow.    In 

particular  the  shed  vorticity  panel  itself  depends  on  yk  to  determine  its  distributed 

vorticity  (yw)k,  this  in  turn  causes  changes  to  (q.)k  and  yk<    Moreover,  the  downstream 

core  vortices  that  constitute  the  wake  are  convected  under  the  influence  of  the  free 

stream  and  the  singularity  distributions  on  the  airfoil  surface  panels  including  the  shed 

vorticity   panel.     The   problem   is    thus    seen    to    be    coupled    from   this    analytical 

standpoint.     Putting   this   in   simple   mathematical   terms,    the    algebriac    system   of 

equations   (Equation   2.18),   representing   the    flow   tangency   conditions    and    Kutta 

condition  for  steadv  flow,  are  no  lonser  linear  because  [he  coefficients  a.,  in  the  left- 

i] 

hand-side  matrix,  are  not  constants  anvmore.  Thev  are  function  of  q.  and  v  instead. 
The  presence  of  non-linearity  is  indeed  what  drives  the  solution  scheme  into  an 
iterative  type  for  unsteady  flow. 


31 


3.  Boundary  Conditions 

We  next  investigate  whether  our  unsteady  flow  model  is  sufficiently 
represented,  before  we  could  proceed  to  search  for  a  numerical  iterative  solution,  by 
matching  the  unknowns  with  the  available  boundary  conditions  at  time-step  L  .  Recall 
that  we  have  introduced  three  more  unknowns  (7w)k,  \  and  0.  in  addition  to 
(q.)k  (j=  1,2,. ...,n)  and  yk-  We  have,  however,  so  far  only  identified  an  extra  boundary 
condition,  namely  the  Helmholtz's  theorem  (Equation  3.2)  in  conjunction  with  the  flow 
tangency  conditions  at  the  n  panel  control  points  and  the  Kutta  condition  of  pressure 
equilibrium  at  the  trailing  edge  panels.  Clearly  we  are  in  deficit  of  two  additional 
conditions  before  attempting  further  endeavour  to  solve  the  entire  system.  Basu  and 
Hancock  [Ref.  3]  suggested  the  following  assumptions  : 

•  The   shed  vorticity  panel  is   oriented  in  the   direction   of  the  local  resultant 
velocity  at  the  panel  mid  point. 

•  The  length  of  the  shed  vorticity  panel  is  proportional  to  the  magnitude  of  the 
resultant  velocity  at  the  panel  mid  point  and  the  step  size  of  the  time-step. 

Following  these  assumptions  thus  permits  us  to  formulate  two  more  boundary 

conditions  as  follows, 


tanOk  =    ,,  w,k  (eqn  3.3) 


\  "  (\  ~  <m)  -J  1  (Uw)k2  +  (VwV  1  («m  3-4) 

where  (Uw)k  and  (Vw)k  are  the  total  velocity  components  in  x  and  y  directions  of  the 
airfoil-fixed  coordinate  system. 

The  flow  tangency  conditions  are  still, 

[(V")^  =  0,   i=l,2,....,n  (eqn  3.5) 

However,  the  Kutta  condition  must  now  include  the  rates  ot  change  or" 
potential  at  the  trailing  edge  panels  (unsteady  Bernoulli's  equation)  which  can  be 
related  directly  to  the  rate  of  change  of  total  circulation.  By  using  a  backward  finite 
difference  approximation  for  this  rate  of  change  of  total  circulation,  we  express  the 
Kutta  condition  as  shown  in  Equation  3.6. 
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KVMk2   "   [(At2  =   2  I  ^  «Pa   "   *1>  it  *    2  <  —  k 

=  2€  Yl<~7k-1  (eqn3.6) 

B.       RIGID  BODY  MOTION  AND  FRAME  OF  REFERENCE 

Consider  a  rigid  airfoil  executing  a  time-dependent  motion,  comprising  linear 
translation  and  angular  rotation  about  the  leading  edge  in  an  inviscid  incompressible 
medium.  We  can  describe  this  arbitrary  motion  at  any  time  instant  tk  as  the  vector  sum 
of  a  mean  velocity  —Vqoi  a  time  dependent  translational  velocity  —  [U(t)  i  4-  V(t)  j] 
and  a  rotational  velocity  ~£2(t)  where  i  &  j  are  unit  vectors  in  the  airfoil-fixed 
coordinate  system  as  shown  in  Figure  3.2. 

If  we  continue,  as  in  steady  flow,  to  determine  the  flow  with  reference  to  the  (x,y) 

coordinate  system  fixed  on  the  airfoil,  an  observer  sitting  on  this  frame  of  reference 

sees  an  unsteadv  stream  velocity,  V  ,        .  made  up  by  the  vector  sum  of  a  mean 

*         stream  r      J 

velocity  Vjq,  a  time  dependent  translational  velocity  [U(t)  i  -  V(t)  j]  and  a  rotational 
velocity  £S(t).  Therefore  in  this  frame  of  reference,  unlike  the  previous  case  where  the 
airfoil  is  allowed  to  move  only  with  constant  linear  velocity,  the  flow  is  still  unsteady  in 
that  Vstream  is  time  dependent. 

Stream  =  Yoo  +  E(U(t)  i  +  V(t)  j)]  +  Q(t)  (y  i  -  x  j)  (eqn  3.7) 

We  redefine  our  disturbance  potential  to  include  the  potential  contributions  (pw 
and  (pcv  from  the  shed  vorticity  panel  and  the  wake  core  vortices  respectively.  Thus  : 

We  then  write  the  total  velocity  with  respect  to  this  frame  of  reference  as, 

Vtotal   "    Stream   +    7(P  fan  ^9) 

Notice  that  this  total  velocity  is  obviously  NOT  the  absolute  velocity  with  respect 
to  an  inertial  coordinate  system.  Such  an  inertial  coordinate  system  will  be  the  one 
where  an  observer  only  sees  an  on-coming  flow  of  V^  with  constant  magnitude  and 
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Figure  3.2     Frame  of  Reference  for  Unsteady  Flow 
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direction.  We  have  to  make  this  distinction  clear  because  in  our  model  on  convection 
of  core  vortices,  we  break  up  the  caiuiation  into  two  steps;  we  first  convect  the  core 
vortices  using  the  resultant  absolute  velocity  with  respect  to  an  inertial  coordinate 
system,  followed  by  determining  their  positions  with  coordinates  relative  to  the  airfoil- 
fixed  axes. 

The  unsteady  flow  Bernoulli's  equation  for  the  pressure  coefficients  on  the  airfoil 
surface  must  be  written  with  respect  to  the  airfoil-fixed  coordinate  system  also.  Giesing 
[Ref.  4]  showed  this  to  be  written,  in  our  notation,  as  : 

S  -  V^Tl  -  (^^")2  "  (^-)2  "  4-1^-  Ceqn3.10) 

p         ^PV00  V00  ^00  voo      At 

where  V  ,        ,  and  V,  ,  ,  are  defined  according  to  Equations  3.7  and  3.9. 

stream'  total  -  ^ 

Equations  3.8,  3.9,  and  3.10  can  be  correlated  to  their  counter-parts  in  steady 
flow,  namely  2.6,  2.8  and  2.9  respectively  with  Vstream  of  Equation  3.7  replacing  the 
Vqo  in  Equation  2.8. 

C.       TIME-DEPENDENT  INFLUENCE  COEFFICIENTS 
1.  Definition  of  Time-Dependent  Influence  Coefficients 

The  influence  coefficients,  An..,  A1..,  Bn..  and  B1..,  involving  the  source  and 
vorticity  distributions  described  in  Section  C  of  Chapter  II  are  still  useful.  These  are 
indeed  time-independent  coefficients  since  they  are  functions  of  geometrical  parameters 
which  are  fixed  in  our  rigid  airfoil.  Additional  influence  coefficients  involving  the  shed 
vorticity  panel  and  the  wake  core  vortices  must  be  defined.  These  coefficients  need  to 
be  computed  in  each  time  step  since  their  positions  vary  relative  to  the  airfoil-fixed 
coordinate  system.  For  that  matter,  as  will  be  made  clear  later,  those  influence 
coefficients  involving  the  shed  vorticity  panel  must  also  be  computed  in  every  iteration 
within  each  time  step  for  the  same  reasoning. 

a.   More  A's  and  B's  Influence  Coefficients 

Following  the  notations  used  previously  in  steady  flow,  we  define,  with  the 

use   of  the  k-subscript   to   denote   time-dependency,  additional  influence  coefficients 

involving  uniformly  distributed  singularities  of  source  and  vorticity.  They  are  the  A's 

and  B's  coefficients  : 

•      (Bnin  +  i)k        :  normal  velocity  component  induced  at   the  Ith  panel  control 
point  by  unit  strength  vorticity  distribution  on  the  shed  vorticity  panel  at  time 

V 
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(Bl.    -i-i\        '•  tangential  velocity  component  induced  at  the  1th  panel  control 
point  by  unit  strength  vorticity  distribution  on  the  shed  vorticity  panel  at  time 

V 

(Axn  +  1  -)k       :  x- velocity  component  induced  at  the  shed  vorticity  panel  mid 
point  by  unit  strength  source  distribution  on  the  j^  panel  at  time  l  . 

(Ay  + ,  .).        :  y-velocity  component  induced  at  the  shed  vorticity  panel  mid 
point  by  unit  strength  source  distnbution  on  the  j1^  panel  at  time  t.  . 

(Bxn  +  1  .)k        :  x-velocity  component  induced  at  the  shed  vorticity  panel  mid 
point  by  unit  strength  vorticity  distribution  on  the  j^  panel  at  time  tk. 

(Byn  +  1  .V        :  y-velocity  component  induced  at  the  shed  vorticity  panel  mid 


point  by  unit  strength  vorticity  distribution  on  the  jm  panel  at  time  t,  . 

(Axh-)k  :  x-velocity  component  induced  at  the  center  of  the  h^  core 

vortex  by  unit  strength  source  distribution  on  the  f*  panel  at  time  L  . 

(Ayh-)k  :  y-velocity  component   induced   at   the   center  of  the  hth   core 

vortex  by  unit  strength  source  distribution  on  the  Ith  panel  at  time  t . 

(B'\.)k  :  x-velocity  component  induced   at   the   center  of  the   h^   core 

vortex  by  unit  strength  vorticity  distnbution  on  the  j1*1  panei  at  time  L- 

(B-v  .).  :  v-velocitv  comDonent   induced   at  the  center  of  the  11th   core 

ill''   K.  •/ 

vortex  by  unit  strength  vorticity  distnbution  on  :he  jril  panei  at  time  L . 

(Bxhn  +  1)k       :  x-velocity  component  induced   at   the  center  of  the   h^  core 
vortex  by  unit  strength  vorticity  distribution  on  the  shed  vorticity  panel  at  time 

V 

(B-vhn+,)k       :  y-velocity  component  induced   at   the  center  of  the   h1*1   core 
vortex  by  unit  strength  vorticity  distribution  on  the  shed  vorticity  panel  at  time 

V 

b.   New  C's  Influence  Coefficients 

The  presence  of  discrete  core  vortices  in  the  wake  requires  the  definition  of 
new  influence  coefficients  involving  point  singularity.  They  are  the  C's  coefficients  in 
our  familiar  notations  : 

•  (Cn;rn)v  :  normal  velocity  component  induced  at   the  Ith  panel  control 
point  by  unit  strength  m     core  vortex  at  time  tk. 

•  (Clim)k  :  tangential  velocity  component  induced  at  the  Ith  panel  control 
point  by  unit  strength  rrr1  core  vortex  at    i:\ 

1      (Cxn^ir         :  x-velocicy  component  induced  at  the  shed  vorticity  panel  mid 
point  by  unit  strength  mth  core  vortex  at  time  l. 

•  (Cyn  +  lnA      :  y-velocity  component  induced  at  the  shed  vorticity  panel  mid 
point  by  unit  strength  nr*1  core  vortex  at  time  t.  . 

•  (Cxhm)k  :  x-velocity  component  induced   at   the  center  of  the   W*1  core 
vortex  by  unit  strength  mth  core  vortex  at  time  L . 
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*      ^h   K  :  y"vei°ci1:y  component   induced   at   the  center  of  the  h^   core 

vortex  by  unit  strength  nr1  core  vortex  at  time  t. . 

2.  Computation  of  Time-Dependent  Influence  Coefficients 

(Bninj_ .).   and  (B^  n  +  1)k  are  computed  exactly  the  same  way  as  Bnr  and  Blv 

are  computed  using  Equations  2.14  and  2.15  with  subscript  n+  1  replacing  j.  Similarly, 

(Ax  -fiOir  and  (Ax..)k  are  calculated  using  Equation  2.12  with  0j  set  to  zero  and 

subsript  i  appropriately  replaced.  Also  (Ayn+lj)k  and  (Ayhj)k  are  calculated  using 

Equation  2.13  with  8.  set  to  zero  and  subsript  i  appropriately  replaced.    We  do  the 

same  for  (Bxn+1  \  and  (Byn^.1  .).   using  Equations  2.14  and  2.15  respectively  and  so 

on  for  all  the  rest  of  A's  and  B's  coefficients.    The  C's  coefficients  will  be  computed 

with  different  expressions  from  those  of  A's  and  B's  because  they  are  the  velocities 

induced  by  unit  strength  core  vortex.  It  can  be  shown  easily,  from  the  geometry  of 

Figure  3.3,  that  their  expressions  take  on  the  foilowmg  forms, 


_    cosje,  -  (9m)k)l 


(C\A  =  "  — -'-,  <ein  3-n> 


2"  (hJ* 


(C'imV  "    -  ,'_,..     .  W  3.12) 


where 


(rim)k  =  VtCxmj-x^  +  Cymj-yJ2] 
xmi  =   1/2(xi  +  xi+1) 
yi^  =  1/2(yi  +  yi  +  1) 


xm  =  x  coordinate  of  m^  core  vortex  at  time  t 


ym  =  y  coordinate  of  mth  core  vortex  at  time  t, 


a  =  arctan(-^±i_Zi_) 


xi  +  l_xi 


vm.  —  v 
(6m\  "  arctanf     '     'm  \ 

xmi~ym 


37 


Bv  the  same  token.  (Cx  _,     V  and  (Cx.     ).    are  comouted  bv  Equation  3.11 

a  —  L  ,m/  k  v       am '  k.  -        * 

while  (Cy  a.  |     ),   and  (Cyhm)k  are  computed  by  Equation  3.12  if  6;  is  set  equal  to  zero 
and  the  subscript  1  appropriately  replaced. 

D.       NUMERICAL  SOLUTION  SCHEME 

1.  The  Flow  Tangency  Conditions 

The  flow  tangency  conditions   of  Equation   3.5   can   be  written   using   the 
influence  coefficients  as  follows, 


n 


I  I  A°ij  (%\  1  +  rk  I  B-y  +■  [  (V^  •  nj  |k 
l-J  .   i=> 

m=  1 

where  (V  ,        ).  is  evaluated  bv  Eauation  3.7  at  the  Ith  panel  control  point 

x     stream' i  -  ^  r 

This  equation,  though  it  seems  complex,  says  nothing  more  than  summing  to 
zero  ail  the  velocity  contributions  due  co  individual  singularity.  Substituting  (Yw)k  from 
Equation  3.2,  collecting  like  terms  and  rearranging  the  equation  into. 

i-1  j-1 

(-[(Vs.eam)iik-ni]-(e/Wk.,(Bni,n+1)k 

-  1 1  <C",A  <rm.!  - Tm)  1 }   .    i"  1.2 n  (eqn  3.14) 

m=  1 

2.  The  Iterative  Solution  Procedure 

Equation    3.14    is    Arranged    in    this    form    because    we    intend    to    solve 

;  ■  ..2 n)  in  "erms  ofyL-  Expanding  Equation  3.14  for  i  =  1,2 n  results  in 

following  matrix  equation, 

[A]{q}w  =  v.    {B}.    +  {  C  }.  (eqn  3.15) 
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Figure  3.3     Influence  Coefficients  due  to  Point  Singularities. 


where  [  A  ]  is  an  n  x  n  matrix  whose  elements  are  known  constants.  {  B  }.   and  {  C  }k 
are  n  x  1  column  vectors  whose  elements  are  known  only  if  the  shed  vorticity  panel  at 
t^  is  established;  that  is,  if  Ak  and  0k  are  known,  then  we  can  calculate  ail  the 
influence  coefficients  on  the  right-hand- side  of  Equation  3.14.  We  therefore  make  use 
of  this  idea  to  formulate  our  iterative  solution  procedure  as  follows  : 

(1)  Project  the  wake  core  vortices  downstream  according  to  the  time  step  and  the 
local  resultant  velocities  at  their  respective  centers  with  respect  to  an  inertial 
coordinate  system. 

(2)  Compute  the  coordinates  of  these  core  vortices  relative  to  the  airfoil-fixed 
coordinate  system  due  to  its  time-dependent  motion. 

(3)  Start  iteration  cycle  for  current  time  step  by  initially  assuming  some  guess 
values  of  Ak  and  0k.  We  can  use,  except  for  the  first  time-step,  values 
obtained  at  previous  time  step.  Compute  then  the  influence  coefficients  needed 
in  Equation  3.14  or  3.15. 

(4)  Obtain  (q.)k  in  terms  of  yk  by  solving  Equation  3.15  as  a  linear  system  with 
two  right-hand-sides  by  the  same  Gaussian  elimination  algorithm  used  in 
steady  flow. 

(5)  Calculate  the  tangential  velocities  at  the  trailing  edge  panels,  all  in  terms  of 

(6)  Invoke  the  Kutta  condition  of  Equation  3.6  (with  some  efforts  in  algebnac 
manipulation)  to  solve  for  yk  since  it  is  the  only  unknown  in  that  equation. 

(7)  Once  yk  is  solved,  (q.)k  are  then  known.  We  can  then  calculate  the  velocity 
components  (Uw)k  and  ( Vw)k  at  the  mid  point  of  the  shed  vorticity  panel. 

(8)  Equation  3.3  and  3.4  hence  enable  us  to  update  the  values  of  Ak  and  ©k. 

(9)  Repeat  the  iteration  cycle  from  steps  (3)  to  (9)  until  converged  values  of  Ak 
and  ©k  are  obtained.  Alternatively  convergence  can  be  set  for  (L"w)k  and 
(Vw)k  instead. 

(10)  Compute  the  tangential  velocities  and  disturbance  potential  at  all  panel 
control  points  in  order  to  determine  the  pressure  distribution  which  can  be 
integrated  to  give  forces  and  moments. 

(11)  Compute  the  resultant  velocities  which  occur  at  the  centers  of  all  the  core 
vortices  that  will  be  convected  down-stream.  These  velocities  must  be  the 
absolute  velocities  with  respect  ro  an  inertial  coordinate  system. 

3.  Computation  or  Velocities 

The    iterative    procedure    mentioned    in    the    previous    subsection    requires 

calculation  of  tangential  velocities  at  the  trailing  edge  panels  and  the  absolute  velocity 

components  (Uw)k  and  (Vw)k.   They  are  computed  differently  due  to  the  use  of  a 

different  frame  of  reference. 
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a.    Tangential  Velocities  on  Airfoil  Panels 

The  tangential  velocities  [(Vl).]k  (i=  1,2 n)  at  all  the  panel  control  points 

are  calculated  using  the  airfoil-fixed  coordinate  frame  of  reference  as  follows  : 


n  n 

'J 


[(vU  ■  I  [  At,i  &h  i  +  \  I B' 
i-i  1-1 


+  [<VstreamVti]k+(Yw)1c(Bt,.n+l)k 

k-1 

+  S  [  (C^k  (rm-l  "  rm)  1    •  1=  ^ >n  ^n  3-16) 

m=l 

b.    Core  Vortices  Convection  Velocities 

The  resultant  velocities  at  all  core  vortices  are  calculated  using  the  inertial 
frame  of  reference  fixed  with  respect  to  Vqq  but  resolving  them  into  components  in  the 
directions  coincide  with  the  airfoil-fixed  coordinate  system  as  shown  below  : 

i-1  1-1 

+  (Voo-i)k+(TA<B\„  +  i)k 

m=l 
m*  h 

n  n 


(Vh)k  -    E  (  (Ayk  (qp,  ]  +  Yk  E  (BV* 

1-1  1-1 

+  <v=o  •i\+(Yw)k(B,h.n+.\ 


k-1 


+  E[(c\m)k(rm.i-rm)i  (•qas.is) 


m=l 
m*  h 
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Notice  the  use  of  V^  instead  of  V  f         in  Equations  3.17  and  3.18.  Also 

-^  stream  n 

the  subscript  h  is  just  an  index  usable  for  any  core  vortex.  We  can  obtain  (U  ).  and 
(V  ).  if  h  is  replaced  by  n-i- 1  in  these  equations. 

4.  Disturbance  Potential  and  Pressure  Distribution 
a.    Why  We  Need  the  Disturbance  Potential 

The  concept  of  disturbance  potential  (p  has  been  instrumental  in  the 
formulation  of  both  the  steady  and  unsteady  flow  problems.  However,  it  has  never 
gone  beyond  using  it  merely  as  a  vehicle  to  understanding  the  superposition  of  simple 
flows.  The  disturbance  potential  need  not  be  solved  for  at  ail  in  the  steady  flow 
problem  formulation.  This  is  because  what  one  really  is  going  after  is  the  spatial 
derivative  of  this  disturbance  potential,  i.e.  the  disturbance  induced  velocity,  from 
which  the  pressure  distribution  can  be  obtained.  We  have,  in  all  our  solutions  so  far, 
been  successful  in  avoiding  any  disturbance  potential  calculation  since  the  concept  of 
influence  coefficients  allows  us  a  direct  evaluation  of  the  velocity.  Unfortunately,  as 
can  be  seen  in  Equation  3.10,  when  we  proceed  further  to  compute  the  pressure 
distribution  on  the  airfoil  surface  in  unsteady  flow,  we  are  faced  with  the  oroblem  of 
evaluating  the  disturbance  potential  <p,  or  more  precisely  the  rate  of  change  of  (p,  which 
we  approximate  by  using  a  backward  finite  difference  expression.  Therefore,  the 
pressure  coefficients  at  the  Vth  panel  control  point  can  be  rewritten,  in  terms  of  non- 
dimensional  variables,  as, 

Lk      Lk-1 

where  (V1).  is  calculated  by  Equation  3.16  and  (V  am)j  is  the  non-dimensional  (by 
Vqq)  form  of  Equation  3.7  evaluated  at  the  i1*1  panel  control  point. 

We  thus  need  to  calculate  at  each  time  step,  the  disturbance  potential  at  all 
the  panel  control  points.  Short  of  having  to  solve  the  Laplace's  equation  by  a  finite 
difference  scheme,  we  evaluate  the  disturbance  potential  (p  by  integrating  rhe  velocity 
field  in  CWO  stages  from  upstream  .it  infinity  to  che  airfoil  leading  edge,  then  along  the 
airfoil  surface  from  the  leading  edge  to  each  panel  control  point.  Care  must  be  taken 
here  to  include  only  the  velocity  contribution  due  to  disturbances. 

One  important  question  arises,  in  this  approach,  as  to  what  value  of 
disturbance  potential  we  should  use  at  infinity  before  we  carry  out  the  line  integral.  We 
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must  therefore  analyse  the  behaviour  of  (p  at  infinity  by  examining  the  singularities  that 
constitute  the  disturbance.  They  are  the  source  and  vorticity  distributions  on  the  airfoil 
surface  and  the  core  vortices  in  the  wake.  These  singularities  induce  no  velocity  at 
infinity  from  the  knowledge  of  simple  flows.  In  other  words,  the  disturbance  potential 
cp  at  infinity  is  independent  of  spatial  coordinates.  The  next  question  we  should  ask  is 
whether  (p  at  infinity  is  time-dependent?  Let  us  adopt  the  view-point  that  if  we  are  at 
infinity  looking  at  our  airfoil  and  its  associated  wake,  we  simply  see  a  point  vortex  with 
a  total  circulation  T0  at  time  tQ.  We  have  already  identified  that  T0  remains  constant 
by  Helmholtz's  theorem.  It  only  gets  redistributed,  as  time  progresses,  over  the  airfoil 
surface  and  in  the  wake.  Notice  that  the  previous  statement  regarding  what  one  would 
see  at  infinity  said  nothing  about  the  source  distributions.  The  source  distributions 
though  vary  (or  get  redistributed  )  as  the  time  progresses,  the  total  source  strength 
necessarily  remains  zero  at  all  time  in  order  to  enforce  a  closed  contour  representing 
the  airfoil  thickness.  This  is  also  the  reason  why  the  unsteady  flow  solution  needs  an 
additional  model  to  handle  the  vorticity  conservation  since  the  source  conservation  is 
already  implicitly  so  for  a  closed  contour  to  exist.  From  these  discussions,  we  are 
certain  that  the  disturbance  potential  (p  at  infinity  is  an  absolute  constant  (independent 
of  time  and  spatial  coordinates)  whose  value  is  fixed  only  by  the  initial  condition  we 
decide  to  start  solving  the  unsteady  problem.  The  actual  value  of  (p  at  infinity  is  in  fact 
immaterial  so  long  as  we  know  it  is  constant  because  its  value  disappears  conveniently 
as  we  subtract  (<Pj)i..i  from  (q>j)k  in  Equation  3.19. 
b.   Computation  of  Disturbance  Potential 

We  begin  by  choosing  an  arbitrary  straight  line  extending  upstream  to 
infinity  from  the  leading  edge  of  the  airfoil  along  a  direction  parallel  to  Vqq.  For 
practical  purposes,  we  set  infinity  at  say  ten  chord  lengths  away  from  the  leading  edge 
since  the  velocities  induced,  at  field  points  thereafter,  by  the  disturbances  are  small 
enough  to  be  negligible.  This  line  is  divided  into  z  panels  with  element  lengths  near  the 
leading  edge  comparable  to  the  panel  sizes  used  on  the  airfoil.  However,  the  panel  size 
is  progressively  increased  to  take  advantage  of  the  inversely  decaying  induced  velocities 
at  larger  distances.  We  compute  the  tangential  components  of  :he  induced  velocities  at 
the  mid  points  of  these  panels  using  influence  coefficients  analogous  to  those  used  on 
the  airfoil  panels.  Using  subscript  f  to  denote  these  panel  mid-points,  we  can  define 
influence  coefficients  (A^,  (Atfjn  +  1)k,  (B1^,  (Blfjn  +  1)k,  and  (Clfm)k  and  compute 
them  using  the  same  expressions  for  calculating  the  A's,  B's  and  C's  coefficients  used 
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before  with  cosG.  replaced  by  ( —  cosa),  sinG.  replaced  by  (  —  since)  and  subscript  i 
replaced  by  f.  With  the  help  of  these  influence  coefficients,  the  tangential  velocity 
induced  by  disturbances  at  the  1th  panel  mid  point  is  : 


I<VVA=    I  I  (AVk  (<&  J  +  Vk  I  (B'fi)k 

1-1  1-1 

+  (Yw)k  (BVn+  ,)k  +  I  [  (C'rm)k  frm.,  -  Tm)  i  (eqn  3.30) 

m=  1 

valid  for  f=  l,2,....,z.  The  disturbance  potential  at  the  airfoil  leading  edge  is  the  sum  of 
the  products  of  the  disturbance  induced  velocity  at  each  panel  and  the  panel  length. 


(<p,e)k  =  ~  I  KvVAc  J  [  (xf  + \ "  xf)2  +  (yf+- 1 "  yf)2  ]  (eqn  3-21) 

f=i 

Similarly,  for  the  line  integral  over  the  airfoil  surface,  we  compute  the 
tangential  component  of  the  disturbance  induced  velocity  at  the  1th  panel  control  point 
using  the  following  equation  : 


KVJk  =    E  t  A«,  (qjV  1  +  Yk  S  B' 
l-i  i-l 


n 


+  <VwV  <*u+ A  +  I !  (cWk  <rm-i  -  rm>  1  <e<in  3-22> 

m=l 

which  is  valid   for  i=  1,2,.... ,n.   Performing   the  line  integration   by  summation,   the 
disturbance  potential  at  the  i     nodal  point  on  the  airfoil  is  : 

i-l 

«Pnode  i>k  ■  CflA  ~  I  :(VVj]k  rj,i  ^  1    '  f0f  n  >  '  *  He 

j_ile 

Me"1 

-  «P,A  "  I  KV^  fjJ  +l   ,        for  i,e  >  i  >   1  (eqn  3.23) 

j-i 

where  r, ,  + .  denotes  the  panel  length. 
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Finally,  the  disturbance  potential  at  the  1th  panel  control  point  is, 

«Pi>k  *    *  I  «Pnode  Ok  +  ^node  t+ 1\  1    '  1=  1'2-'D  ^D  3'24> 

5.  Computation  of  Forces  and  Moments 

The  Cc,  C.  and  C  about  the  leading  edge  are  calculated  in  exactly  the  same 
way  as  it  is  done  for  the  steady  flow  problem  by  integrating  the  pressure  distribution 
(See  section  D-4  of  Chapter  II). 

E.       FLOW  MODELING  OF  SHARP  EDGE  GUST  FIELD 

The  unsteady  flow  solution  described  so  far  can  be  extended  to  the  study  of 
airfoils  penetrating  a  sharp  edge  gust  by  modifying  the  boundary  conditions  with  the 
assumption  that  the  gust  front  remains  straight  while  passing  through  the  airfoil.  The 
same  assumption  has  been  used  in  both  [Ref.  3]  and  [Ref.  4].  An  additional  model  in 
[Ref.  3]  using  distribution  of  singularities  along  the  gust  front  had  successfully 
attempted  to  simulate  the  distortion  of  the  gust  front  passing  over  the  airfoil  surface. 
It  was  shown  that  the  pressure  distributions,  during  the  time  when  the  gust  front 
remained  on  the  airfoil  surface,  were  affected  only  at  the  neighbourhood  of  the  gust 
front.  The  overall  pressure  upstream  and  downstream  of  the  gust  front  stayed 
essentially  the  same.  The  distorted  gust  front  model  is  not  used  in  program  U2DIIF. 
The  use  of  the  relatively  simple  yet  sufficiently  accurate  model  of  a  straight  gust  front 
affords  the  modifications  to  the  unsteady  flow  solution  to  be  confined  only  to  the  flow 
tangency  conditions.  That  is  to  say,  the  expression  of  Vstream  in  Equation  3.7  would 
include  the  gust  velocity  for  panels  that  are  already  in  the  gust  field  during  the 
penetration  phase.  Similarly,  the  computation  of  core  vortex  velocities  using 
Equations  3.17  and  3.18  have  the  gust  velocity  added  to  the  V^q  if  the  core  vortices  are 
already  in  the  gust  field. 

In  an  attempt  to  generalise  the  solution  for  cases  where  airfoils  enter  the  gust 
field  at  an  angle  of  attack,  the  convenient  model  used  in  [Ref.  3]  by  setting  the 
computation  to  proceed,  for  the  undistorted  gust  front  simulation,  so  that  the  gust 
front  always  coincides  with  the  nodal  points  is  difficult  to  implement.  At  any  one  time, 
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if  an  airfoil  enters  a  gust  field  at  an  angle  of  attack.,  the  gust  front  would  appear  in 
between  two  nodes  of  a  particular  panel  on  one  surface  while  the  gust  front  proceeds 
from  node  to  node  on  the  other  surface.  We  therefore  further  modify  the  flow 
tangency  condition  only  on  that  particular  panel  where  the  gust  front  lies  in  between 
two  nodes  by  taking  the  gust  velocity  on  that  panel  to  be  proportional  to  the  fraction 
of  panel  length  partially  submerged  in  the  gust  fieid. 
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IV.  DESCRIPTION  OF  COMPUTER  CODE  U2DIIF 

A.       PROGRAM  U2DIIF  STRUCTURES  AND  CAPABILITIES 
1.  Restrictions  and  Limitations 

The  numerical  formulations  of  both  the  steady  and  unsteady  flow  problems 
outlined  in  the  previous  Chapters  are  coded  in  a  FORTRAN  computer  program  called 
U2DIIF  (See  Appendix  A  for  the  program  listings).  The  present  solution  methods 
treat  the  inviscid  and  incompressible  flow  as  an  approximation  to  the  real  flow  so  long 
as  the  viscous  effect  is  negligible  and  the  flow  stays  attached  on  the  airfoil  surface  at  all 
time.  These  restrictions  are  no  strangers  to  any  one  who  is  familiar  with  any  other 
potential  flow  solution  methods.  Other  than  the  implicit  restrictions  of  potential  flow 
solution,  the  method  is  entirely  general  in  that  the  shape  of  airfoil  is  arbitrary  and  any 
arbitrary  continuous  motion  of  the  airfoil  could  be  simulated  using  either  the  closed 
form  (i.e.  explicit  equations)  or  discrete  data  points  to  describe  the  time-history  of  the 
translational  and  rotational  velocities. 

The  storage  of  the  computer  that  carries  out  the  calculations  may  be  the  other 
limitation  one  should  consider.  The  storage  requirements  grow  rapidly  with  the  number 
of  panels  (n)  and  the  number  of  computation  time  steps  (m).  By  far  the  prime 
contributor  to  this  storage  requirement  comes  from  the  massive  amount  of  influence 
coefficients.  The  number  of  influence  coefficients  increases  with  the  square  of  the 
number  of  panels  (n2).  Each  time  step  increment  adds  (2n  +  m2)  more  influence 
coefficients  due  to  the  formation  of  shed  vorticity.  The  current  program  fixes  the 
maximum  number  of  airfoil  panels  to  200  and  the  maximum  allowable  time  steps  is 
also  200. 

An  additional  constraint  worth  mentioning  concerns  the  gust  field  simulation 
whereby  the  current  solution  methodology  is  valid  except  in  the  use  of  the  same 
pressure  equation  arising  from  the  unsteady  Bernoulli's  equation  (See  Equation  3.10). 
The  fundamental  assumption  underlying  the  derivation  of  this  equation  is  the 
irrotationality  of  the  flow  field.  There  is  no  doubt  that  the  flow  fields  upstream  and 
downstream  of  the  gust  front  are  irrotational.  However,  when  one  needs  to  obtain  the 
pressure  on  the  airfoil  surface,  an  implicit  integration  is  done  across  the  gust  front.  A 
flow  field  inclusive  of  the  gust  front  is  rotational  since  the  line  integral  of  velocity  in  a 
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closed  path  does  not  vanish  when  the  gust  front  is  crossed.  Failing  the  proper 
derivation  of  a  new  pressure  equation  applicable  to  unsteady  rotational  flows,  care 
must  be  exercised  to  regard  the  present  method  as  an  approximate  solution  to  gust 
fields  of  weak  strengths  only. 

2.  Current  Structures  of  U2DIIF  MAIN  Program 

The  overall  program  Logic-flow  chart  is  as  shown  in  Figure  4.1.  The  program 
first  reads  in  the  input  data  from  filecode  1  and  sets  up  the  airfoil  panel  nodes  and 
slopes.  Immediately  after  that,  the  steady  flow  calculations  are  executed  for  the  initial 
angle  of  attack  a.  according  to  the  solution  scheme  described  in  Section  D  of 
Chapter  II.  The  steady  flow  solution  is  included  primarily  to  : 

•  Provide  the  necessary  initial  parameters  for  the  unsteady  flow  solution  to 
proceed  in  time.  In  other  words,  the  steady  flow  solution  handles  the  V^  and 
initial  angle  of  attack  a.  one  decides  to  begin  the  unsteady  flow  calculation. 

•  Allow  the  code  to  function  directly  as  a  steady  [low  solver  as  and  when 
necessary  without  having  to  do  the  time  consuming  unsteady  tlow  iterative 
solution  and  approach  the  steady  flow  as  time  approaches  infinity. 

The  program  terminates  once  the  steady  flow  calculations  are  done  if  the 
program  determines,  based  on  the  input  data  set  by  user,  no  requirement  for  unsteady 
flow  solutions.  Otherwise  the  unsteady  flow  calculations  will  be  activated  by  selecting 
and  computing  the  rigid  body  motions  of  the  airfoil  and  the  corresponding 
computation  time-step  size.  Currently,  all  the  time  dependent  motions  are  equation- 
generated,  they  are  the  positions  and  rates  of  the  transiational  and  rotational  motions. 
Incorporated  as  case-runs  within  the  program  U2DIIF  are  the  following  motions  : 

(1)  Step  change  in  angle  of  attack  from  any  initial  value. 

(2)  Modified- ramp  change  in  angle  of  attack  about  any  pivot  point  from  any 
initial  value. 

(3)  Harmonic  transiational  motion  at  any  angle  of  attack. 

(4)  Harmonic  rotational  motion  about  any  pivot  point  at  any  mean  angle  of 
attack. 

(5)  Sharp  edge  gust  penetration  at  any  angle  of  attack. 

Should  one  decide  to   generate   the  airfoil's  motion  using  discrete  data   points    . 
function  of  time,  the  program  could  be  easily  modified. 

The  computation  time-step  sizes  for  the  harmonic  transiational  and  rotational 
motions  are  constant  values  determined  by  the  frequencies  (FREQ)  and  the  number  of 
computation  per  cycle  (DTS).  For  the  case  of  step  change  in  angle  of  attack,  the 
computation  time-step  size  is  progressively  increasing,  from  a  starting  value  (DTS),  as 
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Figure  4.1     Flow  Chart  for  U2DIIF  Computer  Code. 


49 


Time  Steo  increment 


ves 


Compute  Rigid  Body  Motions  and 
Convect  Core  Vortices  Downstream 


yes 


Comoute  Inr-uence  Coefficients 


Solve  (q-)^  in  terms  of  yk  by 
Gaussian  Elimination  with  2  R.H.S 


Invoke  Kutta  Condition,  Solve  yk 


Update  (U^.)k,    (vw)k,  Ak  and  0k 


Compute  Velocity  [(Vc).],,  Disturbance 
Potential  (<P;)k  and  Pressure  [(Ojl^ 


Compute  Aero  Coeff  (Cg)k,  (Cd)k  and  (Cm)k 


ves 


no  ii 


.cnrou-e   .\esuitant   velocities   or    .ore    /oriicss 


Re-initialising  Parameters   Before 
Starting   the   Next  Time   Step 


Figure  4.1     .  (cont'd.) 


50 


time  increases.  The  case  of  the  modified  ramp  change  in  angle  of  attack,  adopts 
initially  a  constant  computation  time-step  size  (DTS)  during  the  transient  nsing  of  the 
angle  of  attack.  Once  the  final  angle  of  attack  is  reached,  the  computation  time-step 
size  is  progressively  increased  also.  Similarly  for  the  case  of  airfoil  penetrating  a  gust 
field,  the  computation  time-step  size  (DTS)  is  constant  during  the  period  when  the  gust 
front  remains  on  the  airfoil  surface  but  progressively  increases  once  the  entire  airfoil  is 
submerged  in  the  gust  field.  These  variations  in  computation  time  steps  are  to  provide 
greater  flexibility  both  in  capturing  transients  and  covering  relatively  large  total  time  of 
computation  without  having  to  contend  with  the  storage  space  requirements  described 
previously.  These  variations  in  time-step  sizes  described  so  far  are  associated  with 
setting  the  input  parameter  TADJ  to  zero.  If  TADJ  is  chosen  to  be  non-zero,  all  the 
case-runs  would  compute  initially  using  the  starting  time-step  sizes,  based  on  DTS  for 
non-oscillatory  motions  and  FREQ  &  DTS  for  harmonic  motions,  and  the  program 
would  prompt  for  an  user  choice  of  time  step  adjustment.  If  the  answer  is  yes,  the 
program  would  back-track  the  previous  solution  and  recompute  the  current  solution 
using  an  adjusted  time-step  size  that  is  TADJ  times  the  initial  value  i'DTST.  This  special 
time  step  variation  feature  gives  the  program  added  capability  of  allowing  an 
interactive  time  step  selection  during  the  progress  of  unsteady  flow  computation.  The 
ability  to  back-track  and  recompute  the  current  solution  using  a  different  time-step  size 
enhances  the  possibility  of  using  program  U2DIIF  together  with  a  viscous  flow  solver 
forming  an  Inviscid- Viscous-Interactive  solution  scheme  which  often  requires  such  time 
step  variations. 

The  MAIN  program  performs  the  iterative  solution  procedures  set  out  in 
Section  D  of  Chapter  III.  The  convergence  check  during  the  iterative  solution  is  done 
through  the  user  specified  tolerance  between  successive  iterative  solutions  of  both 
(Uw)k  and  (Vw)k.  The  solution  continues  into  the  next  time-step  by  selecting  the  time 
step  size  according  to  the  particular  case-run  and  projecting  all  the  wake  core  vortices 
downstream  so  that  their  new  positions  relative  to  the  airfoil  at  the  new  time  step  can 
be  correctly  determined. 

B.       DESCRIPTION  OF  SUBROUTINES 
1.  Subroutine  BODY 

This  subroutine  is  called  by  subroutine  SETUP  if  the  user  selects  an  airfoil 
that  is  either  a  NACA  XXXX  or  230XX  type.  It  in  turn  calls  subroutine  NACA45  to 
obtain  the  airfoil  thickness  and  camber  distributions  and  returns  with  the  computed 
(x,y)  coordinates  of  the  panel  nodal  points. 

51 


2.  Subroutine  COEF 

This  subroutine  is  called  by  the  MAIN  program  in  the  unsteady  flow 
calculations.  It  utilises,  at  each  iteration  cycle,  the  influence  coefficients  generated  by 
subroutine  INFL  to  calculate  the  coefficients  of  the  matrix  Equation  3.15  by  expanding 
Equation  3.14.  These  matrix  coefficients  are  necessarily  set  up  in  this  way  so  that  the 
source  strengths  could  be  solved  in  terms  of  the  vorticitv  strength  bv  subroutine 
GAUSS  as  a  linear  system  with  two  right-hand-sides. 

3.  Subroutine  COFISH 

This  subroutine  is  called  by  the  MAIN  program  to  set  up  the  coefficients  of 
the  matrix  system  of  Equation  2.18  for  steady  flow  where  the  source  strengths  and 
vorticity  strength  are  solved  simultaneously  by  subroutine  GAUSS  as  a  linear  system 
with  one  right-hand-side.  The  matrix  coefficients  are  calculated  using  Equations  2.16 
and  2.17. 

4.  Subroutine  CORVOR 

This  subroutine  is  called  by  the  MAIN  program  at  nearing  the  end  of  the 
unsteadv  flow  calculations  before  starting  a  new  time  step.  It  comDiites  the  resuitant 
convective  velocities  for  all  the  wake  core  vortices  with  respect  to  an  inertia!  frame  of 
reference  according  to  Equations  3.17  and  3.18  where  all  the  appropriate  influence 
coefficients  are  linked  through  common  block  from  subroutine  INFL. 

5.  Subroutine  FANDM 

This  subroutine  is  used  in  both  the  steady  and  unsteady  flow  calculations.  It 
is  called  by  the  MAIN  program  immediately  after  the  pressure  distribution  over  the 
airfoil  panels  are  known  so  that  it  can  perform  the  simple  integration  of  pressure  in  the 
appropriate  directions  to  give  the  aerodynamic  force  and  moment  coefficients  of  lift, 
drag  and  pitching  moment  about  the  leading  edge  according  to  Equations  2.22  through 
2.26. 

6.  Subroutine  GAUSS 

This  subroutine  is  the  standard  linear  system  solver  that  employs  the  well 
known  Gaussian  elimination  with  partial  pivoting  ind  operates  simuitaneousiv  an  a 
user  specified  number  of  nght-hand-siues.  It  is  called  by  the  MAIN  program  in  ooth 
the  steady  and  unsteady  flow  calculations.  In  order  to  use  GAUSS,  the  coefficients  of 
the  augmented  matrix  must  be  set  up  so  that  GAUSS  will  return  the  solutions 
replacing  the  corresponding  columns  of  the  augmented  matrix  that  were  initially 
occupied  by  the  right-hand-sides.  The  coefficient  set-ups  are  done  by  subroutines 
COFISH  and  COEFF  respectively  for  the  steady  and  unsteady  flow  problems. 
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7.  Subroutine  INDATA 

This  subroutine  is  called  by  the  MAIN  program  to  read  in  the  first  three  sets 
of  data  cards  and  returns  to  the  MAIN  program  if  IFLAG  *  0.  Otherwise  it  continues 
to  read  in  the  fourth  data  card  as  the  NACA  number  corresponding  to  the  type  of 
airfoil  and  calculates  the  thickness  parameters  that  will  be  used  by  subroutine 
NACA45. 

8.  Subroutine  INFL 

This  subroutine  is  the  generator  for  all  the  influence  coefficients  that  need  to 
be  stored  and  used  by  many  subroutines  associated  with  the  unsteady  flow  calculations. 
It  utilises  the  known  relative  geometrical  parameters  of  the  singularities  to  carry  out 
computations  based  on  Equations  2.12  through  2.15,  3.11  and  3.12.  The  MAIN 
program  calls  this  subroutine  in  every  iteration  cycle  of  each  time  step  so  that  the  time- 
dependent  coefficients  can  be  updated  as  and  when  necessary.  Time-independent 
coefficients  are  computed  only  once  in  the  entire  unsteady  flow  solutions.  Those 
influence  coefficients  involving  the  wake  core  vortices  are  updated  in  each  time  step 
while  those  involving  the  shed  vorticity  panel  are  calculated  as  frequently  as  the 
number  of  iterations  take  to  terminate  a  converged  solution.  It,  however,  does  not 
compute  and  store  those  influence  coefficients  needed  for  the  determination  of 
disturbance  potential  (Equation  3.20)  simply  because  they  are  used  only  once  in  each 
time  step. 

9.  Subroutine  KUTTA 

This  subroutine  is  called,  in  the  unsteady  flow  calculations,  by  the  MAIN 
program  during  every  iteration  cycle  in  each  time  step  to  invoke  the  Kutta  condition 
for  unsteady  flow  expressed  in  Equation  3.6.  It  calculates  the  tangential  velocities  at 
the  trailing  edge  panels  using  Equation  3.16  in  terms  of  the  unknown  vorticity  strength 
that  is  manipulated  and  solved  algebraicly. 
10.   Subroutine  NACA45 

This  subroutine  is  called  by  subroutine  BODY  if  the  airfoil  selected  belongs  to 
the  families  of  NACA.  4-digits  airfoils  or  the  NACA  5-digits  airfoils  of  type  230XX 
who  share  common  thickness  distributions  with  the  4-digits  airfoils  having  the  same 
thickness  to  chord  ratio.  The  thickness  and  camber  distribution  data  of  these  airfoils 
are  calculated  and  returned  to  BODY. 
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11.  Subroutine  PRESS 

This  subroutine  is  called  by  the  MAIN  program  to  calculate  the  pressure 
distribution  over  the  airfoil  paneis  after  the  iterative  solution  for  the  unsteady  flow 
problem  has  successfully  met  the  convergence  criterion.  It  first  computes  the 
tangential  velocities  at  all  panel  control  points  using  Equation  3.16,  then  performs  the 
disturbance  potential  evaluation  at  the  current  time  step  according  to  Equations  3.20 
through  3.24.  Together  with  the  disturbance  potential  data  obtained  from  the  previous 
time  step,  it  calculates  the  pressure  distribution  using  Equation  3.19. 

12.  Subroutine  SETUP 

This  subroutine  sets  up  the  panel  nodal  coordinates  for  MAIN  program  by 
reading  the  4th  and  5th  data  sets  of  the  input  file  if  IFLAG  =  1  is  set.  It  skips  the  data 
reading  if  IFLAG  =  0  and  proceeds  to  set  up  the  node  distribution  and  call  subroutine 
BODY  to  calculate  the  airfoil  coordinates.  The  node  distribution  adopts  a  cosine 
formula  in  order  to  have  closelv  oacked  oanels  towards  the  leading  and  trailing  edges 

J  *  i.  W  WW 

for  improvements  in  solution  accuracy.  Regardless  of  how  the  nodal  coordinates  are 
obtained.  SETUP  determines  the  panei  slopes  and  airfoil  perimeter  length. 

13.  Subroutine  TEWAK 

This  subroutine  is  called  by  the  MAIN  program  at  every  iteration  cycle  of 
each  time  step  of  the  unsteady  flow  calculations  to  compute  the  resultant  velocity 
components  at  the  mid  point  of  the  shed  vorticity  panel  using  Equation  3.17  and  3.18 
These  velocity  components  are  necessary  to  ensure  the  correct  establishment  of  the 
shed  vorticity  panel  length  and  orientation  which  governs  the  successful 
implementation  of  the  iterative  solution  scheme  for  the  unsteady  flow  problems. 

14.  Subroutine  VELDIS 

This  subroutine  returns  to  the  MAIN  program  the  velocities  and  pressure 

distributions   for   steady   flow   calculation   using   Equations    2.20   and    2.21.    It   also 

performs  the   evaluation   of  the   disturbance   potential  at   the  panel   control  points. 

Though  these  disturbance  potential  data  are  not  necessary  for  steady  flow  solution, 

iceued  in  the  first  ;;me  step  of  the  unsteady  flow  pressure  calculation. 

C.       INPUT  DATA  FOR  PROGRAM  U2DIIF 

Program  U2DIIF  reads  its  input  data  from  filecode  1.  An  example  of  the  input 
data  file  is  as  shown  in  Appendix  B  for  the  case  where  the  airfoil  nodal  coordinates  are 
input  by  user.  User  could  however  let  the  program  generate  the  nodal  coordinates  if 
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the  airfoil  chosen  happens  to  belong  to  the  family  of  NACA  4-digits  or  5-digits  of  type 
230XX.  To  do  this,  simply  change  I  FLAG  to  zero  in  the  first  item  oi"  the  3rd  set  of 
data  card  and  replace  the  nodal  coordinates  data  in  the  4th  and  3d1  sets  of  data  cards 
by  a  single  data  card  containing  only  the  particular  airfoil's  NACA  number  using 
Format  (15).  Figure  4.2  contains  an  itemised  description  of  the  sequential  input 
variables. 

D.       OUTPUT  DATA  FROM  PROGRAM  U2DIIF 

Appendix  C  contains  a  sample  output  data  generated  by  using  the  input  data  set 
shown  in  Appendix  B.  Due  to  the  repetitive  nature  of  output  as  the  computation  time 
progresses,  only  data  at  selective  time  steps  are  shown.  The  output  data  file  begins  with 
writing  out  what  the  program  has  read  from  the  input  data  file  followed  by  the 
computed  nodal  coordinates  only  if  they  are  program  generated,  otherwise  proceeds  to 
write  the  computed  airfoil  perimeter  length.  The  next  set  of  output  data  are  the  steady 
flow  solution  parameters  of  distributed  source  strengths,  vorticity  strength,  pressure 
and  velocity  distributions  as  well  as  the  force  and  moment  coefficients.  The  output  data 
terminates  at  this  point  unless  unsteady  flow  solution  is  required.  It  then  prints,  for 
each  time  step,  the  unsteady  flow  solution  parameters  similar  to  the  previous  output 
for  steady  flow  with  additional  information  pertaining  to  the  rigid  body  motions  and 
trailing  wake  vortices  data.  An  explanation  of  the  output  variable  names  are  listed  in 
Figure  4.3.   All  output  parameters  are  non-dimensional  quantities. 
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Data  Set  #1 
ITITLE 

Data  Set  #2 
TITLE 

Data  Set  #3 
I  FLAG 

NLOWER 
NUPPER 

Data  Set  #4 
X(I) 

Data  Set  #5 
Y(I) 

Data  Set  #6 
ALPI 

DAL? 

TCON 

FREQ 
PIVOT 

UGUST 

VGUST 

Data  Set  #7 
DELHX 

DELHY 

PHASE 

Data  Set  #8 
TF 
DTS 

TOL 
TAD  J 


Format  (15)  -  I  data  card 

-  Number  of  title  cards  to  be  used  in  Data  Set  42. 

Format  (20A4)  -  ITITLE  data  cards 

-  Headings  to  be  printed  on  output  for  case  run  identification. 

Format  (315)  -  1  data  card 

-  0  if  airfoil  is  NACA  XXXX  or  230XX. 

-  1  otherwise. 

-  Number  of  panels  used  on  airfoil  lower  surface. 

-  Number  of  panels  used  on  airfoil  upoer  surface  (need  not 
be  the  same  as  NLOWER). 

Format  (6F10.6)  if  IFLAG=  1  -  variable  data  cards 

-  x-nodal  coordinates  (divided  by  the  chord  length,  c).  A  total 
of  n+  1  nodal  points  divided  into  6  points  per  data  card. 

Format  (6F10.6)  -  variable  data  cards. 

-  v-nodal  coordinates  ('divided  bv  c)  corresponding  to  the 
Data  Set  #4  if  I  FLAG  =1. 

Format(7F10.6)  -  1  data  card 

-  Initial  angle  of  attack  (AOA)  in  deg. 

-  Increment  in  AOA  in  des  for  non-oscillatory  motions. 

-  Maximum  amoiitude  of  AOA  change  in  deg' for  rotational 
harmonic  motions. 

-  Non-dimensional  rise  time  (V^tjc)  of  AOA  for  motion 
involving  modified-ramp  change  in  AOA. 

-  Non-dimensional  oscillation  frequency  (coc/Vqq)  for 
harmonic  motions. 

-  The  length  from  the  pivot  point  to  the  leading  edge  divided 
by  c  (postive  aft)  for  rotational  motions. 

-  Magnitude  of  gust  velocity  (divided  by  Vqq)  along  Vqq. 

-  Magnitude  of  gust  velocity  (divided  by  Vqq)  perpendicular  to  Vqq 

Format  (3F10.3)  -  1  data  card 

-  Amplitude  of  chordwise  translational  oscillation  divided  by  c. 
(positive  forward). 

-  Amplitude  of  transverse  translational  oscillation  divided  by  c. 
(positive  downward). 

-  Phase  angle  in  de2  between  the  chordwise  and  transverse 
translational  oscillation  with  the  latter  as  reference. 

Format  (4F10.3)  -  1  data  card 

-  Final  non-dimensional  time  to  terminate  unsteady  flow  solution. 

-  Starting  cime-steD  size  for  non-oscillatorv  motiom  DJ  =  " 

-  Number  of  computation  step;    >et  cycle  for  narmomc  motions. 

-  Baseline  time-step  size  for  all  motions  if  TADJ*0.0. 

-  Tolerance  criterion  for  checking  the  convergence  between 
successive  iterations  of  (Uw)k  and  (Vw)k 

-  Factor  by  which  DTS  will  be  adjusted. 


Figure  4.2     List  of  Input  Variables. 


56 


TK 

TKM1 

ALPHA(T) 

OMEGA(T) 

U(T) 

V(T) 

NITR 

VXW 

VYW 

Wake 

THETA. 

GAMMA 

J 

X(J) 

Q(J) 
CP(J) 

V(J) 

CD 
CL 
CM 

M 

X(M) 
Y(M) 
CIRC 


Time  step  tk> 

Time  step  Cj. 

Angle  of  attack  at  time  t. . 

Rotational  velocity  (positive  counter  clockwise)  at  time  tk. 

Chordwise  translation  velocity  (po stive  forward)  at  time  t. . 

Transverse  transiational  velocity  (positive  downward)  at  time  t . 

Iteration  number. 

Iterative  solution  of  (U   ),  . 

Iterative  solution  of(Vw)k. 

Iterative  solution  of  shed  vorticity  panel  length  Ak. 

Iterative  solution  of  shed  vorticitv  oanel  orientation  0,  . 

Iterative  solution  of  the  strength  of  vorticity  distribution. 

Panel  number. 

x-coordinate  of  the  mid  point  of  j     panel. 

Strength  of  source  distribution  on  the  '^  panel. 

Pressure  coefficient  at  the  mid  point  of]™  panel. 

Total  tangential  velocity  at  the  mid  point  of  j     panel 
referenced:  to  the  airfoil-fixed  coordinate  svstem. 

Drag  coefficient. 

Lift  coefficient. 

Pitching  moment  coefficient  about  the  leading  edge. 

Trailing  wake  core  vortex  number. 

x-coordinate  of  the  center  of  m^  core  vortex. 

y-coordinate  of  the  center  of  m^  core  vortex. 

Circulation  strength  of  the  m*  core  vortex. 


Figure  4.3 


List  of  Output  Variables. 
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V.  RESULTS  AND  DISCUSSIOiNS  ON  CASE-RUNS 

This  Chapter  presents  the  results  of  numerous  case-runs  of  U2DIIF  code  for  the 
purpose  of  verifying  the  code.  The  various  airfoils  used  in  the  case-runs  are  deliberately 
chosen  to  be  the  same  as  those  airfoils  where  direct  comparison  of  results  can  be  made 
with  either  theoretical  analyses  and/ or  experimental  data  available  in  the  literature. 

A.       STEP  CHANGE  IN  ANGLE-OF- ATTACK 

1.  Case- Run  Definitions 

Consider  an  airfoil  initially  at  zero  angle  of  attack  to  the  free  stream  Vqq  that 
undergoes  a  step  change  in  angle  of  attack  at  tQ.  The  resulting  flow  should  then 
provide  the  time-dependent  information  on  the  build-up  of  aerodynamic  forces  and 
moments  on  the  airfoil  resembling  the  classical  results  of  Wagner  [Ref.  5]  calculated 
based  on  linearised  theory.  Although  Wagner  prescribed  a  slightly  different  initial 
condition  in  that  the  airfoil  is  initially  at  rest  and  impulsively  started  at  an  angle  of 
attack  and  velocity  Vqq,  the  difference  is  insignificant,  especially  for  a  symmetrical 
airfoil.  This  is  because  the  seemingly  different  initial  conditions  when  translated  into 
the  mathematical  model  means  that  the  step  change  in  AOA  uses  non  zero  initial 
circulation  Tq  at  tQ  with  non  zero  initial  disturbance  potential  at  infinity  if  the  airfoil  is 
cambered.  For  a  symmetric  airfoil,  these  initial  values  are  all  zeroes  and  therefore 
mathematically  would  be  the  same  as  the  initial  conditions  prescribed  by  Wagner. 

2.  Results  and  Discussions 

a.    Von  Mises  8.4%  Thick  Symmetrical  Airfoil 

A  8.4%  thick  symmetrical  Von  Mises  airfoil  is  used  for  this  case-run  where 
the  airfoil  performs  a  0.1  rad  (or  5.73°)  step  change  in  AOA.  Figure  5.1  illustrates  the 
changes  in  the  pressure  distributions  over  the  airfoil  at  time  instances  corresponding  to 
the  airfoil  having  traveled  distances,  in  terms  of  chord  length,  of  0.2.  0.5,  1.0,  2.0  and 
30.  The  associated  trailing  wake  patterns  at  these  time  instances  (less  t»°0)  are  shown 
in  Figure  5.2.  The  time-dependent  build-up  of  aerodynamic  coefficients  of  lift,  drag, 
pitching  moment  and  the  circulation  strength  over  a  computation  period  of  two 
traveled  chord  length  are  shown  in  Figure  5.3.  Notice  that  the  lift,  pitching  moment 
and  circulation  results  are  normalised  by  the  respective  steady  state  values  at  the  same 
AOA.  The  apparently  large  initial  loading  on  the  airfoil  shown  in  Figure  5.3  correlates 
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consistently  with  the  results  of  the  Piston  Theory  of  [Ref.  6]  which  predicts  the  starting 
load  on  an  arbitrary  wing  to  be 

CL  =  4<x/M 

starting 

where  M  is  the  Mach  number.  In  the  case  of  an  incompressible  flow  (M  =  0),  the  initial 
loading  would  be  infinitely  large.  The  same  large  initial  loading  was  obtained  by  Kim 
and  Mook  [Ref.  7]  who  used  continuous  vorticities  as  panel  singularities  instead  of  our 
source  and  vorticity  approach.  Perhaps  what  remains  most  surprising  is  that  the  work 
of  Basu  and  Hancock  [Ref.  3]  did  not  predict  this  initial  loading,  although  they  used 
the  same  singularity  distributions  as  U2DIIF  code.  The  initial  large  loading  in  lift  and 
pitching  moment  decreases  rapidly  over  a  short  time  span,  whereby  the  airfoil  traveled 
approximately  one-tenth  chord  length,  before  rising  in  a  manner  parallel  to  the  Wagner 
Function.  The  drag,  however,  continues  to  decrease  monotonically  after  the  initial 
sharp  fall.  The  circulation  rises,  as  continuous  shedding  of  vorticity  takes  place,  slowly 
from  the  initial  condition  of  zero  to  the  asymptotic  steady  state  vaiue  as  time 
approaches  co.  These  results,  disregarding  the  initial  large  loading  associated  with 
incompressible  flow,  are  in  close  agreement  with  the  results  of  [Refs.  3,4,7]. 
b.    Thickness  Effects  on  the  Wagner  Tunction 

In  order  to  more  closely  correlate  the  results  of  U2DIIF  code  to  the 
theoretical  prediction  of  Wagner,  we  performed  the  step  AOA  change  calculations  for  a 
very  thin  (1%  thickness)  NACA  4-digit  symmetrical  airfoil  which  in  reality  should 
represent  a  flat  plate.  The  results  are  plotted  as  shown  in  Figure  5.4.  Shown  also  on  the 
same  Figure  are  the  results  of  the  8.4%  thick  Von  Mises  airfoil  and  a  25.5%  thick 
symmetric  Joukowski  airfoil.  The  initial  loading  falls  off  less  rapidly  for  the  case  of  the 
simulated  flat  plate  as  compared  to  other  thick  airfoils  but  the  subsequent  rise  in  lift 
follows  very  closely  the  Wagner  Function. 
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Figure  5.1     Pressure  Distributions  at  Various  Time  Instances 

Resulting  from  a  0.1  rad  Step  Change  in  AOA  for  a 

8.4%  Thick  Symmetric  Von  Mises  Airfoil. 
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Figure  5.1     .(cont'd.) 
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Figure  5.3     Time-Dependent  Aerodynamic  Parameters 

Resulting  from  a  0.1  rad  Step  Change  in  AOA  for  a 

8.4%  Thick  Symmetric  Von  Mises  Airfoil. 
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B.       MODIFIED  RAMP  CHAiNGE  IN  ANGLE-OF- ATTACK 
1.  Case-Run  Definitions 

The  case  of  a  step  change  in  AOA  can  be  considered  as  an  usefiii  check,  for 
U2DIIF  code  since  a  handful  of  results  from  other  theoretical  analyses  are  available. 
However,  a  step  change  in  AOA  is  practically  not  realisable  since  all  motions,  short  o{ 
having  infinite  velocities,  take  place  over  a  finite  time  span.  The  step  change  in  AOA 
that  is  practically  possible  is  in  fact  some  form  of  ramp  rise  over  a  short  time  span  with 
large  velocity.  Even  so,  due  to  the  inertia  of  the  airfoil,  an  exact  ramp  rise  in  AOA 
does  not  physically  describe  the  actual  motion  of  the  airfoil  since  finite  time  is  also 
involved  before  the  airfoil  could  build  up  its  ramp  velocity.  Same  argument  holds  at  the 
end  of  the  ramp  rise  before  the  airfoil  could  stop  at  the  final  value  of  AOA.  Therefore 
a  so  called  modified  ramp,  with  some  form  of  rounding  at  the  two  ends  of  a  ramp,  is 
more  likely  to  describe  anything  close  to  what  is  physically  achievable.  The  theoretical 
work  of  Homentcovschi  in  [Ref.  S]  considers  the  case  of  a  flat  plate  that  moves  with 
constant  velocity  and  changes  the  incidence  about  the  mid  chord,  through  a  particular 
ramp  fashion,  described  mathematically  as, 


a(t)  = 


f  0     for  t  <  0, 

6a  (3  -  2t/T)t2/t2     forO  <  t  <x 
5a     for  t  >  t 


where  6a  is  the  magnitude  of  the  AOA  change  and  T  is  the  rise  time  for  the  AOA  to 
reach  its  final  value.  This  particular  function,  plotted  as  shown  in  Figure  5.5,  does  in 
fact  describe  such  a  modified  ramp. 
2.  Results  and  Discussions 
a.    Flat-Plate  Case-Run 

Since  the  results   of  [Ref.  8]   serves  as   another  excellent  source   for  the 
verification  of  U2DIIF  code,  the  obvious  thing  to  do  is  to  use  U2DIIF  to  compute  for 
the  case  of  a  flat  plate,  again  simulated  by  the  1%  thick  NACA-0001  airfoil,  executing 
this  modified  ramp  rise  of  0.1  rad  AOA  over  a  rise  time  of  1.5  chord  length.   Fins 
time  is   chosen  simply  to  facilitate  a  direct  comparison  of  results  to  [Ref.  8]  which  used 
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a  rise  time  of  3  half-chord  lengths.  The  results  of  computation  are  shown  in  Figure  5.6 
and  5.7.  Figure  5.6(a)  takes  a  close  look  at  the  build  up  of  lift  during  the  initial  period 
when  the  airfoil  moves  a  distance  of  six  chord  lengths.  The  lift  initially  rises  to  about 
82%  and  then  decreases  to  about  66%  of  the  steady  state  value  during  the  transient 
rise  time.  Thereafter  it  increases  monotonically  in  a  manner  parallel  to  the  Wagner 
Function.  Figure  5.6(b)  is  a  zoom  view  of  the  rather  slow  convergence  of  lift  to  the 
steady  state  value.  It  takes  the  airfoil  to  cover  a  distance  of  around  50  chord  length 
before  the  lift  builds  up  to  almost  99%  of  the  steady  state  value.  The  same  results  were 
obtained  in  the  theoretical  analysis  of  [Ref.  8].  Figure  5.7  shows  a  collection  of  the 
time-dependent  aerodynamic  parameters  resulting  from  this  particular  case-run. 
b.    Thickness  Effects 

The  same  modified  ramp  function  is  used  on  the  8.4%  thick  Von  Mises 
airfoil.  The  resulting  lift-history  plotted  in  Figure  5.8  shows  a  lower  peak  value  of  lift 
during  the  transient  AOA  rise  as  compared  to  the  case  of  a  Hate  plate  though  a  similar 
trend  of  lift  rise  is  obtained.  Figures  5.9  and  5.10  show  the  results  of  pressure 
distributions  and  trailing  wake  patterns  at  various  time  instances.  One  could  directly 
compare  these  Figures  to  the  corresponding  Figures  arising  from  the  step  AOA  change 
calculations  and  see  the  remarkable  differences  in  transient  characteristics  as  a  result  of 
varying  the  prescribed  motions.  Incidentally,  one  should  realise  that  the  non- 
dimensional  rise  time  of  1.5  chord  length  is  a  deceivingly  large  number.  In  fact,  when 
one  converts  this  to  the  real  time  for  an  airfoil  of  10  ft  chord  length  moving  at  a  low 
Mach  number  of  0.2,  the  rise  time  is  indeed  only  of  the  order  of  0.06  sec.  which  for 
practical  purpose  is  close  enough  to  a  step.  Neverthless,  the  transient  part  of  the  lift 
response  is  entirely  governed  by  how  one  prescribes  the  transient  motion. 
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Figure  5.5     The  Modified  Ramp  AOA  Chance. 
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Figure  5.6    Normalised  Lift  C^C^      Resulting  from  a 

Modified  Ramp  AOA  (6a  =  0.1  raa\  t  =   1.5)  about 

the  Mid  Chord  of  a  NACA-0001  Airfoil. 
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(b)  Normalised  Pitching  Moment  C     C         vs  tV-^.c. 
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Figure  5.7     Time-Dependent  Aerodynamic  Parameters  Resulting  from  a 

Modified  Ramp  AOA  (5a  =  0.1  rad.  t  =   1.5)  about 

the  Mid  Chord  of  a  NACA-0001  airfoil. 
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Figure  5.7     .  (cont'd.) 
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Figure  5.3     Normalised  Lift  C|>  Cg      Resulting  from  a 

Modified  Ramp  AOA  <6a  =  0.1  rod.  t  =   1.5)  about 

the  Mid  Chord  of  a  Miscs  8.4%  Thick  airfoil. 
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Pressure  Distributions  at  Various  Time  Instances  Resulting  from  a 
Modified  Ramp  AOA  (6a  =  0.1  rad,  X  ■   1.5)  about 
the  iMid  Chord  of  a  Mises  8.4%  Thick  airfoil. 
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Figure  5.9     .  (cont'd.) 
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Figure  5.9     .  (cont'd.) 
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Figure  5.9     .  (cont'd.) 
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Figure  5.10    Trailing  Wake  Patterns  at  Various  Time  Instances  Resulting  from  a 

Modified  Ramp  AOA  (6a  =  0.1  rad,  t  =  1.5)  about 

the  Mid  Chord  of  a  Mises  8.4%  Thick  airfoil. 
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C.       TRANSLATIONAL  HARMONIC  OSCILLATION 

1.  Case- Run  Definitions 

Although  the  U2DIIF  code  is  capable  of  computing  unsteady  flow  solution 
for  anv  seneral  translational  harmonic  motion  described  by  a  chordwise  and  a 
transverse  components  bearing  a  given  phase  relationship, 

hy(t)  =  6hy  sin(cot) 
hx(t)  =  h\  sin(oot  +  X) 

where  (0  is  the  oscillation  frequency,  X  is  the  phase  angle  between  the  two  oscillation 
components  and  Sh  &.  Sh  are  the  magnitudes  of  chordwise  and  transverse  oscillations 
respectively.  The  case-run  to  be  presented  in  this  section  selects  the  motion  to  consist 
of  only  the  transverse  component,  i.e.  the  heaving  or  plunging  motion.  A  NACA-0015 
airfoil  is  chosen  for  the  case-run.  The  airfoil  is  initially  at  zero  AOA  with  the 
freestream  Vqq  and  performs  the  plunging  oscillation  at  an  amplitude  of  5h  and  a 
reduced  frequency  of  COc/Vqq 

2.  Results  and  Discussions 

Figures  5.11  and  5.12  present  the  results  of  an  airfoil  executing  a  plunging 
motion  at  an  amplitude  of  0.018c  but  with  two  different  reduced  frequencies  of  4.3  and 
17.0  respectively.  These  numbers  are  chosen  to  coincide  with  those  numbers  used  in 
[Ref.  4].  Excellent  correlations  are  obtained.  Notice  from  these  Figures  that  the 
oscillation  frequency  has  a  great  influence  on  the  magnitudes  of  the  aerodynamic 
parameters  due  to  the  formation  of  significantly  different  trailing  wake  patterns  for  the 
same  oscillation  amplitude.  Also  to  note  is  that  the  width  of  the  resulting  trailing  wake 
is  much  larger  than  the  amplitude  of  the  oscillation,  reinforcing  the  fact  that  the 
unsteady  flow  is  strongly  governed  by  the  shed  vorticity  in  the  trailing  wake.  The  lift 
and  pitching  moment  oscillate  at  the  same  frequency  as  the  airfoil  motion  but  slightly 
out  of  phase,  the  phase  differences  vary  with  the  oscillation  frequency.  The  drag  is 
however  oscillating  at  about  twice  [he  frequency  of  the  airfoil  motion  with  a  negative 
mean  vaiue.  indicating  that  the  plunging  action  indeed  generates  some  propulsive 
thrust.  The  same  conclusion  was  arrived  at  in  the  experimental  work  of  Halfman 
[Ref.  9]  using  a  symmetrical  NACA-0012  airfoil. 
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Figure  5.11     Harmonic  Pluncinsz  Motion  of  a  NACA-0015  Airfoil 
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Figure  5.12     .  (cont'd.) 
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D.       ROTATIONAL  HARMONIC  OSCILLATION 

1.  Case- Run  Definitions 

The  case  of  an  airfoil  oscillating  harmonically  in  pitch  will  be  uniquely  defined 
only  if  a  pivot  point  is  prescribed.  If  a  fixed  pivot  point  is  used  in  the  calculations,  take 
for  example  :he  leading  edge,  any  pitching  motion  about  a  pivot  point  other  than  the 
leading  edge  would  need  r.o  be  described  as  composed  of  a  pitching  and  a  translation 
motions  about  the  leading  edge.  Program  U2DIIF  handles  such  conversion 
automatically  without  having  the  user  to  figure  out  the  combined  motion.  This  applies 
aiso  to  the  case  of  modified  ramp  rise  in  AOA.  The  harmonic  pitching  oscillation  is 
described  by, 

a(t)  =  5a  sin(oot) 

where  6a  and  (0  are  the  amplitude  and  frequency  of  oscillation  respectively. 

2.  Results  and  Discussions 

""he  results  for  the  case  of  the  3:4%  thick.  Von  Mises  symmetric  airfoil, 
oscillating  at  an  amplitude  of  0.01  rad  (or  0.573°)  at  a  rather  high  reduced  frequency 
of  (Oc.'V-jq  =  20.0  about  the  leading  edge,  are  shown  in  Figure  5.13.  The  lift,  drag  and 
pitching  moment  time-history  as  well  as  the  trailing  wake  patterns  are  very  much 
similar  to  the  case  of  a  plunging  airfoil  at  frequency  of  the  same  order  of  magnitude. 
The  differences  are  clearly  in  the  magnitudes  and  phase  angles.  These  results  check 
closely  to  those  of  [Ref.  3].  Figure  5.14  shows  the  results  of  the  same  Mises  airfoil 
performing  another  harmonic  oscillation  at  a  lower  reduced  frequency  of  0.8  and  a 
large  amplitude  of  0.3973  rad  about  a  pivot  point  0.5  chord  length  ahead  of  the  leading 
edge.  [Ref.  4]  conducted  the  same  analysis  for  this  pitch  oscillation  although  the 
reason  for  using  such  a  high  amplitude  of  almost  23°  was  not  clear.  It  is  envisaged 
that  such  high  amplitude  may  result  in  flow  separation.  Nevertheless,  the  case-run  is 
carried  out  assuming  validity  of  attached  flow  for  the  sake  of  comparing  the  results. 
Perhaps  in  inherent  advantage,  in  :he  light  of  U2DIIF  -ode  verification,  with  the  use 
jf  implitude  .n  .".his  case-run  is  that  a  discrepancy,   if  any,  would  show  up 

significantly.  Due  to  the  use  of  different  computation  time-step  size  the  pressure 
distributions  on  the  airfoil,  shown  in  Figure  5.14,  do  not  correspond  one-to-one  at 
exactly  the  same  angular  positions  as  those  presented  in  [Ref.  4].  However,  the  angular 
positions  are  matched  to  within  0.001°,  0.05°  and  0.8°  respectively  for  the  three 
pressure  distributions  shown.  They  correlate  very  well  to  the  results  of  [Ref.  4]. 
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(c)  Trailing  Wake  Pattern  at  the  end  of  2  cycles. 
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(d)  Trailing  Wake  Pattern  at  the  end  of  4  cycles. 


Figure  5.13     .  (cont'd.) 
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Figure  5.14     Harmonic  Pitching  Motion  about  a  Pivot  0.5c  ahead  of 

the  Leading  Edge  of  a  8.4%  Thick  Von  Mises  Airfoil 

5a  =  0.3973  rad,  coc/V^  =  0.8. 
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0.3  0j6  0.7  08 

x.'c  (Chordwise  Station  as  Fraction  of  Chord) 


a9 


to 


(b)  Pressure  Distribution  at  a  =    -0.0775  rad. 


Figure  5.14     .  (cont'd.) 
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(c)  Pressure  Distribution  at  a  =  0.1876  rad. 


Figure  5.14     .(cont'd.) 
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E.       SHARP  EDGE  GUST  FIELD  PENETRATION 

1.  Case-Run  Definitions 

The  case  of  an  airfoil  penetrating  a  sharp  edge  gust  fieid  can  be  computed 
using  the  U2DIIF  code  by  specifying  the  components  of  gust  velocity  along  and 
perpendicular  to  the  freestream  Vx  and  the  angle  of  attack  q[  the  airfoil.  The  results 
that  are  presented  here  consider  the  case  of  airfoil  penetrating  a  sharp  edge  vertical 
gust  at  zero  AOA,  in  view  of  generating  information  on  the  time-dependent  lift 
resembling  the  classical  results  of  Kussner  [Ref.  10]  based  on  linearised  theory.  The 
gust  front  is  taken  to  be  at  the  leading  edge  at  tQ  with  only  the  transverse  (vertical) 
component  of  0.25VQQ. 

2.  Results  and  Discussions 

Figures  5.15  and  5.16  demonstrate  the  variation  of  pressure  distributions  and 
trailing  wake  patterns  respectively  during  and  shortly  after  the  gust  front  moves  past 
the  airfoil.  It  is  interesting  to  see  that  the  resulting  wake  pattern  after  the  entire  airfoil 
is  submerged  in  the  gust  field  is  as  if  being  split  by  the  gust  front  into  two  portions 
curling  in  opposite  directions.  [Ref.  3]  predicted  a  similar  behaviour  for  the  case  o:  an 
undeformed  gust  front.  Due  to  the  use  of  the  modified  flow  tansencv  condition  to 
handle  the  situation  when  the  gust  front  happens  to  fall  in  between  two  nodal  points, 
the  pressure  distributions,  predicted  by  U2DIIF  code,  lie  in  between  the  results  of  the 
undeformed  and  deformed  gust  front  models  used  in  [Ref.  3].  We  therefore  conclude 
that  this  modified  flow  tangency  condition  produces  sufficiently  accurate  results 
without  adding  complication  in  deformed  gust  field  modeling  which  in  turns  limits  the 
application  to  only  sharp  edge  gusts.  The  present  method  would  therefore  preserve  the 
generality  for  extending  calculations  to  other  types  of  gust  front.  Another  comparison 
is  made,  as  shown  in  Figure  5.17,  by  plotting  the  build-up  of  lift  as  a  function  of 
distance  traveled  by  the  airfoil  in  chord  lengths.  Shown  in  the  same  Figure  are  the 
Kussner  Function  and  the  results  obtained  from  another  case-run  using  a  25.5%  thick 
Joukowski  airfoil. 
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x/c  (Chordwise  Station  as  Fractic 


(a)  tV^/c  =  0.25 


Figure  5.15     Pressure  Distributions  at  Various  Time  Instances  Resulting 

from  a  8.4%  Thick  Von  Mises  Airfoil  Penetrating  a 

Vertical  Sharp  Edge  Gust  of  0.25  V^. 
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Figure  5.16    Trailing  Wake  Patterns  at  Various  Time  Instances  Resulting 

from  a  S.4%  Thick  Von  Mises  Airfoil  Penetrating  a 
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Figure  5.16    .(cont'd.) 
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Figure  5.17     Normalised  Time-Dependent  Lift  Cf  Cj      due  to 

Airfoils  of  Various  Thicknesses  Penetrating 

a  Vertical  Sharp  Fdge  Gust  of  0.25V  x. 
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VI.  CONCLUDING  REMARKS 

A.  GENERAL  COMMENTS 

The  U2DIIF  computer  code  has  been  developed  for  the  purpose  of 
demonstrating  the  successful  extension  of  the  well  known  panel  methods,  which  have 
been  used  extensively  for  steady  flow  problems,  into  a  powerful  numerical  tool  for 
solving  the  unsteady  flow  problems.  The  mathematical  modeling  of  the  various  types 
oi^  unsteady  flows  has  been  done  with  the  goal  of  preserving  the  generality  of  the 
methods.  The  intention  is  to  present  a  method  that  has  the  minimum  inherent 
limitations  and  restrictions  so  that  its  usage  for  future  applications  and  developments 
could  remain  appealing. 

The  validation  of  the  U2DI1F  code  has  been  done  through  the  various  case-runs 
of  the  numerous  types  of  unsteady  flow  problems.  The  results  of  each  case-run  has 
been  shown  to  be  weil  correlated  to  the  results  obtainable  from  the  literature  in  the 
form  of  theoretical  analyses,  numerical  calculations  based  on  different  variants  of  panel 
singularities  and  in  cases  where  limited  experimental  data  are  available.  The  ability  of 
those  case-runs  using  an  airfoil  as  thin  as  1%  thickness  to  produce  results  that 
correlate  accurately  to  the  theoretical  flat  plate  results  is  perhaps  a  remarkable 
robustness  possessed  by  the  present  unsteady  flow  solution  methods. 

B.  ENHANCING  U2DIIF  PROGRAM'S  CAPABILITY 

It  has  been  noted  in  Chapter  IV  that  the  current  U2DIIF  code  limits  the  total 
number  of  panels  to  200  and  the  total  computation  time  steps  to  200  also.  These  limits 
are  not  at  all  rigid  and  can  be  easily  increased  if  the  computer  storage  space  is  not 
critical.  A  point  to  note  is  that  the  computation  time  will  grow  rapidly  as  these  limits 
are  raised.  The  current  linear  system  solver,  the  Gaussian  elimination  algorithm,  used 
in  the  code  must  be  concurrently  improved  upon  to  efficiently  reduce  the  computation 
time  required  for  the  iterations  in  each  time  step.  A  close  examination  of  the  matrix 
Equation  3.15,  where  the  linear  system  solver  is  needed  in  every  iteration,  reveals  that 
the  coefficients  of  the  left-hand-side  matrix  [  A  ]  are  time-independent  constants. 
Therefore  the  Gaussian  elimination  algorithm  need  only  be  done  once  for  the  entire 
unsteady  flow  calculations  as  far  as  the  left-hand-side  matrix  is  concerned.  One  could 
then   perform,    for   each    iteration,    the    manipulation    of  the    two    right-hand-sides 
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according  to  the  steps  taken  for  the  reduction  of  the  left-hand-side.  This  should  cut 
down  the  computation  time  significantly  against  the  current  method  of  manipulating 
both  the  left-  and  right- hand- sides  simultaneously  in  each  iteration.  The  savings  in 
computation  time  was  not  pursued  in  the  development  of  U2DIIF  code  because  the 
prime  concern  was  to  demonstrate  that  the  basic  iterative  solution  scheme  works  for 
the  unsteady  flow  problems. 

Another  improvement  that  one  may  consider  is  to  enable  the  code  to  be 
continuable  from  a  time  step  where  previous  calculations  were  terminated.  One  sees 
this  requirement  necessary  not  only  m  the  case  of  premature  termination  of 
computation  due  to  some  unforeseen  circumstances,  but  also  if  one  needs  to  prolong 
the  computation  time.  Certainly  with  the  current  code  structure,  one  has  no  choice  but 
to  perform  the  calculation  right  from  the  beginning. 

A  farther  extension  of  U2DIIF  code  to  the  computation  of  unsteady  flow 
problems  involving  multiple  airfoils  may  be  worth  pursuing.  Other  research  works  that 
could  be  done  based  on  U2DIIF  code  are  in  the  area  of  incorporating  more  variety  of 
rigid  boay  motions  into  the  code  either  in  the  form  of  closed-form  equations  or 
tabulated  time  history  of  the  positions  and  rates  of  motions.  It  is  important  that  one 
should  use  as  close  as  possible,  in  the  code,  a  rigid  body  motion  that  describes  the 
physical  motion  before  generating  any  numerical  unsteady  flow  results  for  meaningful 
comparison  to  test  data.  This  fact  has  been  well  illustrated  and  emphasized  in  the 
comparison  of  results  of  case-runs  involving  a  step  change  to  that  of  a  ramp  change  in 
AOA  with  a  fast  rate  which  one  could  regard  as  a  step  in  reality.  However,  remarkable 
difference  in  transient  aerodynamics  has  been  shown. 
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APPENDIX  A 
U2DIIF  PROGRAM  LISTINGS 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCGCCCCCCCCCC 

c  c 

C     PROGRAM  U2DIIF  C 

c  c 

C  UNSTEADY  MOTION  OF  A  TWO-DIMENSIONAL  AIRFOIL  C 

C  IN  INCOMPRESSIBLE  INVISCID  FLOW  C 

C  USING  PANEL  METHODS  BASED  ON  THE  HESS  &  SMITH  C 

C  C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

COMMON  /BOD/  IFLAG, NLOWER,NUPPER,NODTOT,X(202) ,Y(202) , 
+  COSTHE(201) ,SINTHE(201) ,SS,NP1,NP2 

COMMON  /WAR/  VYW , VXW , WAKE , DT 

COMMON  /WAK2/  VYWK,VXWK 

COMMON  /SING/  Q ( 200 ) .GAMMA, QK( 200) ,GAMK 

COMMON  /CORV/  CV(200)  ,XC(200)  ,  YC(200)  ,M  ,TD  ,  CWX(200)  ,  CWY(200 ) 

COMMON  /POT/  PHI (200) ,PHIK(200) 

COMMON  /GUST/  UG{200) , VG(200) ,XGF, UGUST , VGUST 

COMMON  /EXTV/  UE(200) 

DIMENSION  XXC(200),YYC(200) 
C 

C        INPUT  FROM  FILE  CODE  5  AND  SET  UP  PANEL  NODES  AND  SLOPES 
C 

PI      =  3.1415926535 

WRITE  (6,1003) 
1003  FORMAT  (/////,'  DATA  READ  FROM  FILE  CODE  V  ,11) 

CALL  INDATA 

CALL  SETUP 

READ  (1,501)  ALPI, DALP,TCON,FREQ, PIVOT, UGUST, VGUST 

WRITE  (6,501)  ALPI,DALP,TCON,FREQ, PIVOT, UGUST, VGUST 
501   FORMAT  (7F10.6) 

READ  (1,501)  DELHX , DELHY , PHASE 

WRITE  (6,501)  DELHX, DELHY, PHASE 

READ  (1,501)  TF,DTS,TOL,TADJ 

WRITE  (6,501)  TF.DTS,TOL,TADJ 

IF  (IFLAG  .EQ.  0)  WRITE  (6,1005) 
1005  FORMAT  (///  '  COORDINATES  OF  AIRFOIL  NODES', 
+  //,3X,'  X/C* ,6X  '  Y/C ,/) 

IF  (IFLAG  .EQ.  0)  WRITE  (6,1010)  (X(I) , Y(I) , 1=1 ,NODTOT+l) 
1010  FORMAT  (F10.6,F10.6) 

WRITE  (6,1020)  SS 
1020  FORMAT(//,'  AIRFOIL  PERIMETER  LENGTH  =  ' ,F10.6,/) 
C 

C     STEADY  FLOW  CALCULATION  AT  ALPI 
C 

ALPHA   =  ALPI 

WRITE  (6,1030)  ALPHA 
1030  FORMAT  (//,'  STEADY  FLOW  SOLUTION  AT  ALPHA  =  ' ,F10.6,/) 

IF  (ALPHA  .GT.  90.)     GO  TO  200 

COSALF  -   COS(ALPHA*PI/180.) 

SINALF  =  SIN(ALPHA*PI/180.) 

CALL  COFISH(SINALF, COSALF) 

CALL  GAUSS (1) 

CALL  VELDIS( SINALF, COSALF) 

CALL  FANDM( SINALF, COSALF) 
C 

C        INITIALISATION  FOR  UNSTEADY  FLOW  CALCULATION  TO  BEGIN 
C 


HX 

=   0.0 

HY 

=   0.0 

HXO 

=   0.0 

HYO 

=   0.0 
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100 


c 
c 
c 


3 
2 

1 

60 

1051 

40 


C 

c 
c 


51 
50 


: 
c 
c 


c 

c 
c 


33 


DHX 
DHY 
UX 
UY 

ALP  = 
DA 
COSDA 
SINDA 
OMEGA 
XGF 
ANGLE 
COSANG 
SINANG 
DO  100 
UG(IG) 
VG(IG) 
PHA 
VXW 
VYW 
GAMK 
T 
M 
TOLD 


=  0.0 


=  0.0 
=  0.0 
=  0.0 
ALP  I 
=  0 


0 

=  1.0 
=  0.0 
=  0.0 
=  0.0 

=  ALPI*PI/130. 
=  COS (ANGLE) 
=  SIN(ANGLE) 
IG  =  1,N0DT0T 
=  0.0 
=  0.0 
=  PHASE*PI/180. 
=  COSALF 
=  SINALF 
=  GAMMA 
=  0.0 
=  0.0 
=  0.0 


+  ATAN(VGUST/(1.+UGUST)) 


RIGID  BODY  MOTIONS  OF  AIRFOIL 


IF  (FREQ 

IF  (DALP 

IF  (TCON 

ALPHA 

COSALF 

SINALF 
DT      =  DTS 

TD      =  DTS 

GO  TO  60 
IF  ((UGUST  .EQ 
=  DTS 
=  DTS 


NE.  0.0)  GO  TO  1 
EO.  0.0)  GO  TO  2 
NE.  0.0)  GO  TO  3 
ALPI  +  DALP 
C0S(AL?HA*?I/130. 
SIN(AL?HA*PI/130. 


0.0)  .AND.  (VGUST  .EQ.  0.0))  GO 

DT 

TD 

GO  TO  60 
DT 

TD 
T 

WRITE 

FORMAT 
+  i  *** 

M       =  M  +  1 

IF  (T  .GT.  TF)  GO  TO  200 


TO  200 


=  2.0*PI/(FREQ*DTS) 

=  DT 
=  DT 
(6,1051) 

(/////   I   ****************************************   / 

BEGIN  UNSTEADY  FLOW  SOLUTION   ****',/, 


STORE  CORE  VORTEX  COORDINATES  FOR  TIME  STEP  ADJUSTMENTS 


IF 
DO 


g 


XXC ( I ) 
YYC(I)   = 
IF  (FREO 
I?  fDALP 
[F   TCON 


•EQ.  1)  GO  TO 
I  =  1,M-1 


50 


XC(I) 
YC(I) 
NE.  0.0)  GO  TO  11 

,io.  d.o)  30  ro  :: 

.NE.  D.O)  JO  TO 


70 


STEP  CHANGE  IN  AOA 

IF  (TADJ  .NE.  0.0)  GO  TO 
TD  =  FLOAT (M+1)*DTS 
GO  TO  70 


MODIFIED  RAMP  CHANGE  IN  AOA 


IF  (T  .GT.  TCON)  GO  TO  34 
DAL     =  DALP  *  (3.-2.*T/TCON)*(T/TCON)**2 
ALPHA    =  ALPI  +  DAL 
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34 


C 

G 
C 


22 


111 


120 


121 
110 


C 

c 
c 


11 


c 
c 


COSALF 
SINALF 

DA 

COSDA 
SINDA 
OMEGA 
DHX 
DKY 
UY 

MTCON 
GO  TO 
DAL 
ALPHA 
COSALF 
SINALF 
DA 

COSDA 
SINDA 
OMEGA 
DHX 
DHY 
UY 

IF  (TAD J 
TD 
GO  TO  70 


) 


COS(ALPHA*PI/180 

3 IN (ALPHA*?!/ 130 
ALPHA  -  ALP 
COS(DA*?I/130.) 
SIN(DA*PI/130. ) 

-  (DAL?*?I/130.)  *  (6.*T/(TC0N*TC0N))  *  (l.-T/TCON) 
PIVOT  *  (1. -COSDA) 

-  PIVOT  *  SINDA 
PIVOT  *  OMEGA 

H 


70 
=  0.0 
ALP  I 


+  DALP 


COS(ALPHA*PI/180.) 
SIN(ALPHA*PI/180.) 

0.0 
1.0 
0.0 
0.0 
0.0 
0.0 
0.0 

.NE.  0.0)  GO  TO  70 
FLOAT (M+l -MTCON) *DTS 


SHARP  EDGE  GUST  (UGUST  AND/ OR  VGUST) 


XGF 
DO 


110 


=  T 


UG(IG) 


VG(IG) 

XG 

XGP1 

IF  (IG 

IF  (XGF 

IF  (XGF 

FAC 

UG(IG 

VG(IG 

GO  TO 
UG(IG) 

VG(IG) 

GO  TO 
IF  (XGF 

IF  (XGF 

FAC 

UG(IG 

VG(IG 

GO  TO  11 
UG(IG)   = 

VG(IG) 
CONTINUE 

IF  (XGF 

IF  (TADJ 

IF  (XGF 

GO  TO  70 


Y(IG)*SINALF 
Y(IG+I)*SINALF 


11 


11 


IG  =  l,NODTOT 

=  0.0 

=  0.0 

=  X(IG)*COSALF  -f 

=  X(IG+l)*COSAL5 

LT.  NLOWER+1)  GO  TO  120 

.LE.  XG)  GO  TO  110 

.GE.  XGP1)  GO  TO  111 

=  (XGF  -  XG)/(XGP1  -  XG) 

=  UGUST*FAC 

=  VGUST*FAC 

0 

UGUST 
■  VGUST 
0 

LE.  XGP1)  GO  TO  110 
.GE.  XG)  GO  TO  121 
=  (XGF  -  XGP1)/(XG  -  XGP1) 
=  UGUST*FAC 
=  VGUST*FAC 
0 

UGUST 
=  VGUST 


LE.  COSALF)  MGUST  =  M 
.NE.  0.0)  GO  TO  70 
,GT.  COSALF)  TD  =  FLOAT (M+l -MGUST )*DTS 


TRANSLATION  HARMONIC  OSCILLATION 


))  GO  TO  12 

■  SIN(FREQ*T  +  PHA) 


IF  (DALP  .NE.  0.0' 

HX  =  DELHX 

HY  =  DELHY  *  SIN(FREQ*T) 

DHX  =  HX  -  HXO 

DHY  =  HY  -  HYO 

UX  =  DELHX*FREQ*COS(FREQ*T+PHA) 

UY  =  DELHY*FREQ*COS(FREQ*T) 
GO  TO  70 

ROTATIONAL  HARMONIC  OSCILLATION 
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c 

12    DAL  =  DAL?*SIN(FREO*T) 

ALPHA  =  ALPI  +  DAL 

COSALF  =  C0S(ALPHA*PI/18Q. 

SINAL?  =  SIN(AL?KA*PI/I30. 

DA  =  ALPHA  -  ALP 

COSDA  =  COS(DA*PI/130. 

SIMDA  =  SIN(DA*PI/130. 

OMEGA  =  -  (DALP*PI/180'.)  *  FREQ  *  COS(FREQ*T) 

UY  =  PIVOT  *  OMEGA 

DHX  =  PIVOT  *  (1. -COSDA) 

DHY  =  -  PIVOT  *  SIMDA 
C 
C        TRANSFORM  CORE  VORTEX  COORDINATES  W.  R.  T.  NEW  AIRFOIL  POSITION 


C 


70    IF  (M  .SO.  1)  GO  TO  30 

DO  90    I  =  1,M-1 

XC(I)    =  XXC(I)  +  CWX(I)  *  DT 

YC(I)    =  YYC(I)  +  CWY(I)  *  DT 

XCO     =  XC(I) 

YCO     =  YC(I) 

XC(I)    =  XCO*COSDA  -  YCO*SINDA  +  DHX 
90    YC(I)    =  XCO*SINDA  +  YCO*COSDA  +  DHY 
30    CONTINUE 

"•/RITE  (6,1001)  T,DT 
1001  FORMAT  (/////,'  TIME  STEP  TK  =  ' , F10 . 6 , 10X, ' TK  -  TKM1  =  " ,F10.6,/) 


FORMAT  {/////,'    TIME  STEP  TK  =  '  ,F 

WRITS  (6,1004)  ALPHA , OMEGA , UX , UY 
FORMAT  (/,'  ALPHA (T)  =  ',F10.6,5X, 


1004  FORMAT  (/,'  ALPHA (T)  =  ',F10.6,5X,'  OMEGA(T)  =  '  F10.6,/. 

+  '  U(T)      =  ^FIO^SX,'  V(T)      =  ',F10.6,///, 

+  IX,1  NITR    VXW       VYW      WAKE      THETA     GAMMA',/) 
C 

C      CALCULATE  THE  TRAILING  EDGE  WAKE  ELEMENT 
C 

MITR    =  ] 
10    WAKE     =  SQRT(VYW*VYW+VXW*VXW)*DT 
THEMP1   =  ATAN2(VYW,VXW) 
COSTHE(NPl)  =  C0S(THENP1' 
SINTHE(NPl)  =  SIN(THENP1 

WRITE  (6,1002)  NITR, VXW, VYW, WAKE, THENP1,GAMK 
1002  FORMAT  (15 , 4F10 . 6 , E14 . 6) 

X(NP2)  =  X(NP1)  +  WAKE*C0STHE(NP1) 
Y(NP2)  =  Y(NPl)  +  WAKE*SINTHE(NP1) 
CALL  INFL  (NITR) 

CALL  COEF  (SINALF, COSALF, OMEGA, UX,UY) 
CALL  GAUSS (2) 

CALL  KUTTA  (ALPHA , SINALF , COSALF , OMEGA ,UX,UY) 
CALL  TEWAK  (SINALF , COSALF) 
TOL1     =  ABS(VYW  -  VYWK) 
TOL2     =  ABS(VXW  -  VXWK) 

IF  ((TOL1  .LT.  TOL)  .AND.  (TOL2  .LT.  TOL))  GO  TO  20 
VYW     =  VYWK 
VXW     =  VXWK 
IF  (NITR  .GT.  200)  STOP 
NITR    =  NITR  +  1 
GO  TO  10 
20    WRITS  f 6, 1011)  NITR 

1911  FORMAT        ZONVERGED  SOLUTION  OBTAINED  AFTER  NITR  =  ',13) 
:\LL    PRESS  (SINALF  :3SALF,  OMEGA,  UX.'JY) 

IF     JUST  .EQ.  ).0)  .AND.  ' VGUST  .EQ,  0.0))    GO  TO  300 
CALL  FANDM  (.SiNANG,  COSANG) 
GO  TO  400 
300   CALL  FANDM  (SINALF, COSALF) 
400   CONTINUE 
C 

C     ADJUST  TIME  STEP  (TADJ  .NE.  0.0)  IF  NECESSARY 
C 

IF  (TADJ  .EQ.  0.0)  GO  TO  95 
WRITE  (5,2001) 
2001  FORMAT  (//,'  DO  YOU  WANT  TO  ADJUST  TIME  STEP  ?  0  -  NO ,  1  -  YES') 
READ  (5,*)  IDT 
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(IDT  .EQ.  0)  GO  TO  95 
=  TADJ  *  DT 
=  TOLD  +  DT 
'6,1006) 
1006  FORMAT  (//,'  BACK-TRACK  COMPUTATION  AND  ADJUST  TIME-STEP',//) 
10 
C 

C        WAKE  ELEMENT  LEAVES  TRAILING  EDGE  AS  A  CORE -VORTEX 
C 
95    CV(M)    =  SS*( GAMMA -GAMK) 

XC(M)    =  X(NP1)  +  0.5*WAKE*COSTHE(NP1) 
YC(M)    =  Y(NP1)  +  0.5*WAKE*SINTHE(NP1) 
CWX(M)  =  VXW 
CVVY(M)  =  VYW 
WRITE  (6,1052) 
1052  FORMAT  (//,'  TRAILING  VORTICES  DATA',//, 
+  4X, 'M' ,4X, 'X(M) ' ,6X, 'Y(M) ' ,6X, 'CIRC ,/) 
DO  900   I  =  1,M 
900   WRITE  (6,1050)  I,XC(I) ,YC(I) ,CV(I) 
1050  F0RMAT(I5,3F10.6) 

CALL  CORVOR  (SINALF, COSALF) 
C 

C        RE- INITIALISE  PARAMETERS  FOR  NEXT  TIME  STEP  CALCULATION 
C 

DO  30    I  =  l,NODTOT 
0(1)    =  QK(I) 
PHI(I)   =  PHIK(I) 
30    CONTINUE 

GAMMA   =  GAMK 
ALP     =  ALPHA 
HXO     =  HX 
HYO     =  HY 
TOLD    =  T 
DT      =  TD 
T       =  T  +  TD 
GO  TO  40 
200   STOP 
END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C     SUBROUTINE  BODY(Z,SIGN,X,Y)  C 

C  C 

C  RETURN  COORDINATES  OF  POINT  ON  THE  BODY  SURFACE  C 

C  C 

C  Z  =  NODE-SPACING  PARAMETER  C 

C  X,Y  =  CARTESIAN  COORDINATES  C 

C  SIGN  =  +1.  FOR  UPPER  SURFACE  C 

C  -1.  FOR  LOWER  SURFACE  C 

C  C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  BODY(Z , SIGN,X, Y) 
COMMON  /PAR/  NACA,TAU,EPSMAX,PTMAX 
IF  (SIGN  .LT.  0.0)     Z  =  1.  -  Z 
CALL  NACA45 (Z , THICK , CAMBER , BETA) 
X       =  Z  -  SIGN*THICK*SIN(BETA) 
Y       =  CAMBER  +  SIGN*THICK*COS(BETA) 
RETURN 
END 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

C     SUBROUTINE  COEF  (SINALF , COSALF, OMEGA ,UX,UY)  C 


c 


c 


C  SET  COEFFICIENTS  OF  N  EOUS  ARISING  FROM  FLOW  C 

C  TANGENCY  CONDITIONS  AT  MID  POINTS  OF  PANELS  C 

C  SOLVING  THE  N-SOURCE  STRENGTHS  IN  TERMS  OF  THE  C 

C  VORTICITY  STRENGTH  (RESULTING  IN  2  RHS)  C 

C  KUTTA  CONDITION  IS  SATISFIED  SEPARATELY  TO  OBTAIN  C 

C  THE  VORTICITY  STRENGTH  C 

C  THIS  SOLUTION  METHOD  IS  DESIRED  FOR  UNSTEADY  FLOW  C 


C 


C 


109 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  COEF  (SINALF, COS ALF , OMEGA, UK , UY) 
COMMON  /BOD/  IFLAG,NLOWER ,NU??ER, NODTOT , X( 202) , Y(202 ) , 
+  C0STHE(201) /SINTHE(201) ,SS,NP1,NP2 

COMMON  /COF/  A(201,211),NEQS 
COMMON  /SING/  Q (200) , GAMMA, OK (200) ,GAMK 
COMMON  /WAR/  VYW , VXW , WAKE , DT 

COMMON  /CORV/  CV(200) ,XC(200) , YC(200 ) ,M,TD , CCVX(200 ) , CCVY(200) 
COMMON  /INF1/  AAN(201,201) , BBN( 201 , 201 ) ,AYNP1(201) ,3YNP1(201) 
COMMON  /INF2/  CCN(201 , 200) , CCT(201 , 200 ) , CYNP1 (200) , CXNP1 ;200 ) 
COMMON  /GUST/  UG(200) , VG(20Q ) ,XGF , UGUST , VGUST 
NEQS     =  NODTOT 
NP1      =  NODTOT  +  1 
NP2     =  NODTOT  +  2 


C 

c 
c 


c 
c 
c 


90 


110 


c 
c 

c 


c 
c 
c 


INITIALISE  COEFFICIENTS 

DO  90    I  =  1, NODTOT 
DO  90    J  =  1,NP2 
A(I,J)   =  0.0 

SET  LHS  MATRIX  A(I,J) 


DO  120 
XMID 
YMID 
B 

DO  110 
A(I,J) 
B 
CONTINUE 


=  1, NODTOT 

0.5  *  (X(I)  +  X 

0.5  *  (Y(I)  +  Y 

0.0 

=  1, NODTOT 

AAN(I,J) 

B  +  BBN(I/J) 


III)) 


FILL  IN  THE  RIGHT  HAND  SIDE 


A(I,NP1)  =  -B  + 
A(I,NP2)  =  -BBNi 
+  +  SINTHE(I)*  (( 
+  -  COSTHE(I)*  (( 


BBN(I,NP1)*SS/WAKE 

I,NP1)*GAMMA*SS/WAKE 

1 . +UG ( I ) ) *COSALF-VG ( I ) *SINALF+UX ) 

1 . +UG ( I ) ) *S INALF+VG ( I ) *COSALF+UY ) 


+  OMEGA*(YMID*SINTHE(I)  +  XMID^COSTHE (I ) ) 


ADD  CORE  VORTEX  CONTRIBUTION 


IF  (M  .EQ.  1)  GO  TO  140 

MM1      =  M  -  1 

DO  100   N  =  1,MM1 

A(I,NP2)  =  A(I,NP2)  -  CCN(I,N)*CV(N) 
100   CONTINUE 
140   CONTINUE 
120   CONTINUE 

RETURN 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c 
c 
c 
c 
c 
: 
c 
c 
c 
c 


SUBROUTINE  COFISH(SINALF , COSALF) 

SET  COEFFICIENTS  OF  LINEAR  SYSTEM  -  M+l  EQUATIONS 

M  EOUS  -  FLOW  TANGENCY  AT  MID  POINTS  OF  PANELS 
L  EQU   -  ECATTA  :ONDITION  AT  TRAILING  ZDGE  PANELS 

THIS  SOLUTION  METHOD  IS  EFFECTIVE  FOR  STEADY  FLOW.  MO 
ITERATION  IS  REQUIRED,  N-SOURCE  STRENGTHS  AND  1 
VORTICITY  STRENGTH  ARE  SOLVED  SIMULTANEOUSLY 


C 

c 
c 

c 
c 

: 
c 
c 
c 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  COFISH ( SINALF , COSALF) 

COMMON  /BOD/  IFLAG ,NLOWER,NUPPER , NODTOT , X(202) ,Y(202 ) , 
+  COSTHE(201) ,SINTHE(201) ,SS,NP1,NP2 

COMMON  /COF/  A(201,211),KUTTA 
COMMON  /MUM/  PI,PI2INV 
KUTTA    =  NODTOT  +  1 
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FLOG 

FTAN 


C  INITIALISE  COEFFICIENTS 

C 

DO  90    J  =  1,KUTTA 
90    A(KUTTA,J)    =0.0 
C 

C  SET  VN  =  0  AT  MID-POINT  OF  I-TH  PANEL 

C 

DO  120   I  =  l,NODTOT 

XMID    =  .5*(X(I)  +  X(I+1)) 

YMID    =  .5*(Y(I)  +  Y(l+1)) 

A(I/KUTTA)  =0.0 
C 

C  --    FIND  CONTRIBUTION  OF  J-TH  PANEL 

C 

DO  110   J  =  l,NODTOT 

FLOG    =0.0 

FTAN    =  PI 

IF  (J  .EO.  I)     GO  TO  100 

DXJ     =  XMID  -  X(J) 

DXJP    =  XMID  -  X(J+1) 

DYJ     =  YMID  -  Y(J) 

DYJP    =  YMID  -  Y(J+1) 

=  .5*ALOG(  (DXJ?*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ)) 
=  ATAN2(DYJP*DXJ-DXJ?*0YJ,DXJP*DXJ+DYJP*DYJ) 
100   CTIMTJ   =  COSTHE(I)*COSTHE(J)  +  5INTHE (I ) *SINTHE ( J) 

STIMTJ   =  SINTHE{I)*COSTHE(J)  -  COSTHE(I)*SINTHE( J) 

A(I,J)   =  PI2INV*(FTAN*CTIMTJ  +  FLOG*STIMTJ) 

3       =  PI2INV*(FLOG*CTIMTJ  -  FTAN*STIMTJ) 

A(I,XUTTA)  =  A(I,KUTTA)  +  B 

IF  ((I  .GT.  1)  .AND.  (I  .LT.  NODTOT))GO  TO  110 
C 

C  --    IF  I-TH  PANEL  TOUCHES  TRAILING  EDGE, 

C  ADD  CONTRIBUTION  TO  XUTTA  CONDITION 

C 

A(KUTTA,J)  =  A(KUTTA,J)  -  3 

A(KUTTA,KUTTA)  =  A(KUTTA,KUTTA)  +  A(I,J) 
110   CONTINUE 
C 

C  FILL  IN  KNOWN  SIDES 

C 

A(I,KUTTA+1)  =  SINTHE(I)*COSALF  -  COSTHE(I )*SINALF 
120   CONTINUE 

A(KUTTA,KUTTA+1)  =  -  (COSTHE(l)  +  COSTHE(NODTOT) )*COSALF 
+  -  (SINTHE(l)  +  SINTHE(NODTOT))*SINALF 

RETURN 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  SUBROUTINE  CORVOR  (SINALF , COSALF)                               C 

C  c 

C  COMPUTE  THE  LOCAL  VELOCITIES  OF  CORE  VORTICES            C 

C  C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  CORVOR  (SINALF , COSALF) 

COMMON  /BOD/  IFLAG,NLOWER ,NUPPER,NODTOT ,X(202) , Y(202) , 
+  COSTHE(201) ,SINTHE(201) ,SS,NP1,NP2 

COMMON  / S ING /  0(200), GAMMA , OK ( 2 0  0 ) , GAMK 

COMMON  /WAK/  VYW , VXW , WAKE , DT 

COMMON  /'CORV/  CV(200)  rXC(200)  ,  YCC200  )  ,M,TD  ,  CCVX(200  )  ,  CCVY(200  ) 

COMMON  /INF3/  AMY( 200 , 201 ) , BMY( 200 , 201 ) 

COMMON  /INF4/  CMY(200 , 200 ) , CMX(200 , 200) 

COMMON  /POT/  PHI (200) ,PHIK(200) 

COMMON  /GUST/  UG ( 200 ) , VG ( 200 ) , XGF , UGUST , VGUST 

IF  (M  .EQ.  1)  GO  TO  40 

MM1     =  M  -  1 
C 

C        VELOCITY  COMPONENTS  OF  CORE  VORTICES  AT  CURRENT  TIME  STEP 
C 

UGC     =0.0 

VGC     =0.0 
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DO  10   N  =  1,MM1 

KG      =  XC(N)*COSALF  +  YC(N)*SINALF 
I?  (XG  .GT.  XGF)  GO  TO  5 
UGC     =  UGUST 
7GC     ■  7GUST 
5     VY      =  SS*3MY(N,NP1)*(GAMMA-GAMK)/WAKE+ 
+  (l.+UGC)*SINALF+VGC*COSALF 

VX      =  SS*AMY(N/NP1)*(GAMMA-GAMX)/WAKE+ 
+  (l.+UGC)*COSALF-VGC*SINAL 

DO  20    J  =  i,NODTOT 

VY      =  VY  +  AMY  N,J)*OK(J)  +  BMY(N, J)*GAMK 
VX      =  VX  -  3MY(N,J)*6K(J)  +  AMY(N,J)*GAMK 
20    CONTINUE 
C 

C        ADD  CORE  VORTEX  CONTRIBUTION 
C 

DO  30   MC  =  1,MM1 
IF  (MC  .EQ.  N)  GO  TO  30 
VY      =  VY  +  CMY(N,MC)*CV(MC) 
VX      =  VX  +  CMX(N,MC)*CV(MC) 
30    CONTINUE 
C 

C        COORDINATES  OF  CORE  VORTICES  AT  NEXT  TIME  STEP 
C 

CCVX(N)  =  VX 
CCVY(N)  =  VY 
10    CONTINUE 
40    CONTINUE 
50    CONTINUE 
RETURN 
END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c  c 

C      SUBROUTINE  FANDM(SINALF , COSALF)  C 

c  c 

C  COMPUTE  AND  PRINT  OUT  CD , CL , CM  C 

C  INTEGRATE  PRESSURE  DISTRIBUTION  BY  TRAPEZOIDAL  RULE      C 

C  c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  FANDM(SINALF , COSALF) 

COMMON  /BCD/  IFLAG ,NLOWER ,NUPPER,NODTOT ,X(202) , Y(202) , 
+  COSTHE(201) ,SINTHE(201) ,SS,NP1,NP2 

COMMON  /CPD/  CP(200) 

CFX  =0.0 

CFY  =0.0 

CM  =0.0 
DO  100   I  =  l,NODTOT 

XMID  =  .5*(X(I)  +  X(I+1)) 

YMID  =  .5*(Y(I)  +  Y(I+1)) 

DX  =  X(I+1)  -  X(I 

DY  =  Y(I+1)  -  Y(I 

CFX  =  CFX  +  CP(I)^DY 

CFY  =  CFY  -  CP(I)*DX 

CM  =  CM  +  CP(I)*(DX*XMID  +  DY*YMID) 
100   CONTINUE 

CD      =  CFX*COSALF  +  CFY*SINALF 

:l    =  :fy*cosalf  -  cfx*sinalf 

write  6,1000)  :d  :l  ::: 

1000  format  (  /  /  ,  '   cd  =  '  ,  f10  . 6  ,  '    cl  =  '  ,  f10  .6 ,  '   cm  =  '  ,  f10  . 6  ) 

RETURN 
END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C     SUBROUTINE  GAUSS (NRHS)  C 

C  C 

C  SOLUTION  OF  LINEAR  ALGEBRAIC  SYSTEM  BY  C 

C  GAUSS  ELIMINATION  WITH  PARTIAL  PIVOTING  C 

C  C 

C  A  COEFFICIENT  MATRIX  C 

C  MEQNS      =  NUMBER  OF  EQUATIONS  C 
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C  NRHS      =  NUMBER  OF  RIGHT  HAND  SIDES  C 

C  c 

C  RIGHT-HAND  SIDES  AND  SOLUTIONS  STORED  IN  C 

C  COLUMNS  NEQNS+1  THRU  NEONS +NRHS  OF   A  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

SUBROUTINE  GAUSS (NRHS) 

COMMON  /COF/  A(201, 211), NEQNS 

NP      =  NEONS  +  1 

NTOT    =  NEONS  +  NRHS 
C 

C  GAUSS  REDUCTION 

C 

DO  150   I  =  2,NEQNS 
C 

C  —   SEARCH  FOR  LARGEST  ENTRY  IN  (I-l)TH  COLUMN 

C  ON  OR  BELOW  MAIN  DIAGONAL 

C 

IM      =1-1 

I MAX    =  IM 

AMAX    =  ABS(A(IM,IM)) 

DO  110   J  =  I, NEONS 

IF  (AMAX  .GE.  ABS(A(J,IM)))  GO  TO  110 

IMAX    =  J 

AMAX    =  ABS(A(J,IM)) 
110   CONTINUE 
C 

C  --    SWITCH  (I-l)TH  AND  IMAXTH  EQUATIONS 

C 

IF  (IMAX  .EQ.  IM)    GO  TO  140 

DO  130   J  -   IM,NTOT 

TEMP     =  A(IM,J) 

A(IM,J)  =  A(IMAX,J) 

A(IMAX,J)  =  TEMP 
130   CONTINUE 
C 

C  ELIMINATE  (I-l)TH  UNKNOWN  FROM 

C  ITH  THRU  (NEQNS)TH  EQUATIONS 

C 
140   DO  150   J  =  I, NEONS 

R       =  A(J,IM)/A(IM,IM) 

DO  150   K  =  I, NTOT 
150   A(J,K)   =  A(J,K)  -  R*A(IM,K) 
C 

C  BACK  SUBSTITUTION 

C 

DO  220   K  =  NP,NTOT 

A(NEQNS,K)  =  A(NEQNS,K)/A(NEQNS, NEQNS) 

DO  210   L  =  2,NEQNS 

I       =  NEQNS  +  1  -  L 

IP      =1+1 

DO  200   J  =  IP. NEQNS 
200   A(I,K)   =  A(I,K>  -  A(I.J)*A(J,K) 
210   A(I,K)   =  A(I,K)/A(I,I) 
220   CONTINUE 

RETURN 

END 

C  C 

C     SUBROUTINE  INDATA  C 

c  c 

C  SET  PARAMETERS  OF  BODY  SHAPE  C 

C  FLOW  SITUATION,  AND  NODE  DISTRIBUTION  C 

C  C 

C  USER  MUST  INPUT  C 

C  NLOWER  =  NUMBER  OF  NODES  ON  LOWER  SURFACE  C 

C  NUPPER  =  NUMBER  OF  NODES  ON  UPPER  SURFACE  C 

C  PLUS  DATA  ON  BODY  AND  SUBROUTINE  BODY  C 

C  C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
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10 
501 
502 
503 


SUBROUTINE  INDATA 

DIMENSION  TITLE (20) 

COMMON  /BOD/  IFLAG,NLOWER,NUPPER ,NODTOT ,X(202) , Y(202^ 
+  COSTHE(201) ,SINTKE(201) ,SS,NP1,NP2 

COMMON  /PAR/  NACA, TAU, EPSMAX, PTMAX 

READ  (1,501)  ITITLE 

WRITE  (6,501)  ITITLE 

DO  10    I  =  1, ITITLE 

READ  (1,502)  TITLE 
WRITE  (6,503)  TITLE 
FORMAT (31 5) 
FORMAT ( 20 A4) 
FORMAT ( IX, 20A4) 

READ  (1,501)  IFLAG,NLOWER,NUPPER 

WRITE  (6,501)  IFLAG,NLOWER,MUPPER 

IF  (I FLAG  .NE.  0)  RETURN 

READ  (1,501)  NACA 

WRITE  (6,501)  NACA 


IEPS 

IPTMAX 

ITAU 

EPSMAX 

PTMAX 

TAU 

IF  ( 

PTMAX 

EPSMAX 

RETURN 


PS 


NACA/1000 
=  NACA/100  -  10*IEPS 
=  NACA  -  1000*IEPS  - 
=  IEPS*0.01 
=  IPTMAX*0.1 
=  ITAU*0.01 

.LT.  10)  RETURN 
=  0.2025' 
=  2.6595*PTMAX**3 


100*IPTMAX 


END 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  : 

C      SUBROUTINE  INFL  (MITR)  C 

C  C 

C  CALCULATE  INFLUENCE  COEFFICIENTS  C 

c  c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  INFL  (MITR) 

COMMON  /BOD/  IFLAG,NLOWER,NUPPER,NODTOT ,X(202 ) , Y(202) , 
+  C0STHE(201) ,SINTHE(201) ,SS,NP1,NP2 

/NUM/  PI,PI2INV 
/WAK/  VYW,VXW,WAKE,DT 

/CORV/  CV(200) ,XC(200) ,YC(200) ,M,TD , CCVX(200) , CCVY(200) 
AAN(201,201),BBN(201,201),AYNP1(201),BYNP1(201' 
/CCT(201,200),CYNP1(200),CXNP1(200 
,BMY(200,201 
,CMX(200,200' 
•  GT.  0))  GO  TO  510 


C 
C 
C 
C 
C 

c 
c 


COMMON 
COMMON 
COMMON 
COMMON 
COMMON 
COMMON 
COMMON 
IF  ((M 


/INF1/ 
/INF2/ 
/INF3/ 
/INF4/ 
.GT.  1) 


CCN(201,200 
AMY(200,201 
CMY(200,200 
.OR.  (NITR' 


AAN(I,J)  :  NORMAL  VELOCITY  INDUCED  AT  MID-POINT  OF  I-TH  PANEL 

BY  UNIT  STRENGTH  DISTRIBUTED  SOURCE  ON  THE  J-TH  PANEL 

BBN(I,J)  :  NORMAL  VELOCITY  INDUCED  AT  MID-POINT  OF  I-TH  PANEL 

BY  UNIT  STRENGTH  DISTRIBUTED  VORTEX  ON  THE  J-TH  PANEL 


=  l,NODTOT 
.5*(X(I)  + 
, 5*(Y( I j  - 
=  l.NODTOT 


X(I+1) 

Y  I 


+1)) 
H)J 


.EQ 


100 


DO  120 

XMID 

YMID 

DO  . .  3 
.  LOG 
FTAN 
IF  (J 
DXJ 
DXJP 
DYJ 
DYJP 

FLOG    =  .5*ALOG((DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ)) 
FTAN    =  ATAN2(DYJP*DXJ-DXJP*DYJ,DXJP*DXJ+DYJP*DYJ) 
CTIMTJ   =  COSTHE(I)*COSTHE(J)  +  SINTHE ( I )*SINTHE ( J) 
STIMTJ   =  SINTHE(I)*COSTHE(J)  -  COSTHE ( I ) *SINTHE ( J) 
AAN(I,J)  =  PI2INV*(FTAN*CTIMTJ  +  FLOG*STIMTJ) 


0.0 
PI 

■  I) 
XMID 
XMID 
YMID 
YMID 


GO  TO  100 
X(J) 
X 

Y 
Y 


J+l) 

j+i) 
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110 
120 
510 


.35 


C 

c 

c 
c 
c 

u 

c 
c 
c 


130 


c 
c 
c 
c 
c 
c 
c 


140 


BBN(I,J)  =  PI2INV*(FLQG*CTIMTJ  -  FTAN*STIMTJ) 

CONTINUE 
CONTINUE 
CONTINUE 

T 

XMID   =  .5*(x(r)  +  X 

YMID    =  .5*(Y(I)  +  Y 

DO  130 
FLOG 
FTAN 
IF  (J 
DXJ 
DXJP 
DYJ 
DYJP 
FLOG 
FTAN 
CTIMTJ 
STIMTJ 


NP1 

.5*(X(I 
.5*(Y(I 
=  1,NP1 

0.0 

PT 


SO 


XMID 
XMID 
YMID 
YMID 

.5*ai 


i+i 

T+1 


135 


GO  TO 

-  X(J) 

-  X(J+1) 

"  Y  J)  X 

-  Y  J+l) 

OG((DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ)) 

=  ATAN2(DYJP*DXJ-DXJP*DYJ,DXJP*DXJ+DYJP*DYJ) 
=  COSTHE(I)*COSTHE(J)  +  SINTHE (I)*SINTHE ( J) 
=  SINTHE(I)*COSTHE(Jv 


COSTHE(I)*SINTHE(J) 
AAN(I,J)  =  ?I2INV*'(FTAN*CTIMTJ  +  FLOG*STIMTJ) 
BBN(I,J)  =  PI2INV*(FLOG*CTIMTJ  -  FTAN*STIMTJ) 


AYNPl(J) 


BYNPl(J) 


AYNPKJ 

3YNP1(J 
CONTINUE 

DO  140 

XMID 

YMID 

J 

DXJ 

DXJP 

DYJ 

DYJP 

FLOG 

FTAN 

CTIMTJ 

STIMTJ 

AAN(I,J 

BBN(I,J 
CONTINUE 

CCN(I,J) 
CCT(I,J) 


IF  (M  .EQ 

mmi    = 

IF  (NITR 


Y  -  VELOCITY  INDUCED  AT  MID  POINT  OF  WAKE  ELEMENT 
(NP1-TH  PANEL)  BY  UNIT  STRENGTH  DISTRIBUTED  SOURCE 
ON  J-TH  PANEL 

Y  -  VELOCITY  INDUCED  AT  MID  POINT  OF  WAKE  ELEMENT 
(NP1-TH  PANEL)  3Y  UNIT  STRENGTH  DISTRIBUTED  VORTEX 
ON  J-TH  PANEL 


=  PI2INV*(FTAN*COSTHE(J 

=  PI2INV*(FLOG*COSTHE(J 


-  FLOG*SINTHE(J)) 
+  FTAN*SINTHE(J)) 


I  =  l,NODTOT 


X(I 

Y(I 


1  +  1 
1  +  1 


=  .5" 
=  .5* 

••   NP1 

■   XMID 

=  XMID 

-   YMID 

=  YMID 

=  . 5*ALOG ( ( DXJP*DX JP+DYJP*DYJP ) / (DXJ*DX J+DYJ*DYJ ) ) 

=  ATAN2(DYJP*DXJ-DXJP*DYJ,DXJP*DXJ+DYJP*DYJ) 

=  COSTHE(I)*COSTHE(J)  +  SINTHE (I )*SINTHE ( J* 

=  SINTHE(I)*COSTHE(J)  -  COSTHE (I )*SINTHE( j' 

=  PI2INV*(FTAN*CTIMTJ  +  FLOG*STIMTJ) 

=  PI2INV*(FLOG*CTIMTJ  -  FTAN*STIMTJ) 


DO  220 

XMID 

YMID 

DO  210 

DX 

DY 

DIST 

COSTHN 

SINTHN 

CTIMTN 

STIMTN 

CCN(I,N) 


T 


N 


NORMAL  VELOCITY  INDUCED  AT  MID-POINT  OF  I-TH  PANEL 
BY  UNIT  STRENGTH  N-TH  CORE  VORTEX 

TANGENTIAL  VELOCITY  INDUCED  AT  MID-POINT  OF  I-TH  PANEL 
BY  UNIT  STRENGTH  N-TH  CORE  VORTEX 

1)  RETURN 
M  -  1 

GT.  0)  GO  TO  520 
=  1 . NODTOT 
0.5*(X(I)  +  X(I+1 
0.5*(Y(I)  +  Y(I+1 
=  1,MM1 

XMID  -  xc(n; 

YMID  -  YC(N] 

SQRT(DX*DX+DY*DY) 

DX/DIST 

DY/DIST 

COSTHE (I )*COSTHN  +  SINTHE (I )*SINTHN 

SIMTHE(I)*COSTHN  -  COSTHE (I )*SINTHN 

■  -PI2INV*CTIMTN/DIST 
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210 
220 
520 


C 
C 
C 
C 
C 
C 
C 


c 

c 
c 

c 

c 
c 


c 
c 
c 
c 
c 
c 
c 


CCT(I,N) 
CONTINUE 
CONTINUE 
CONTINUE 

XMID 

YNID 

DO  230 

DX 

DY 

DIST 

COSTHN 

SINTKN 

CTIMTN 

STIHTN 

CCN(I,N) 

CCT(I,N) 

CYNPl(N) 
CXNPl(N) 


=  -PI2INV*STIMTN/DIST 


N 


SH 


NP1 

0.5*(X(I)  +  X 

0.5*(Y(I)  +  Y 

=  1,MM1 

XMID  -  XC(N) 

YMID  -  'IC(U) 

SORT(DX*DX+DY* 

DX/DIST 

DY/DIST 

COSTHE(I)*COSTHN  + 

SINTKE(I)*COSTHN  - 

■  -?I2INV*CTIMTN/DIST 

■  -PI2INV*STIMTN/DIST 


fDY) 


SINTHE(I)*SINTHN 
COSTHE(I)*SINTHN 


Y  -  VELOCITY  INDUCED  AT  MID  POINT  OF  WAKE  ELEMENT 
'VNP1-TH  PANEL)  BY  UNIT  STRENGTH  N-TH  CORE  VORTEX 

X  -  VELOCITY  INDUCED  AT  MID  POINT  OF  WAKE  ELEMENT 
(NP1-TH  PANEL)  BY  UNIT  STRENGTH  N-TH  CORE  VORTEX 


CYNPl(N)  =  -PI2INV*C0STHN/DIST 
CXNPl(N)  =  +?I2INV*SINTHN/DIST 


50   CONTINUE 
AMY(N,J) 


3MY(N,J) 


IF  (NITR 

DO  320   N 

XMID     = 

YMID 

DO  310   J 

DXJ 

DXJ? 

DYJ 

DYJP 

FLOG    = 

FTAN 

AMY(N,J' 

BMY(N,J 
310  CONTINUE 
320  CONTINUE 
530   CONTINUE 

DO  330   N 

XMID 

YMID 

J       = 

DXJ 

DXJP 

DYJ 
;? 

FLOG    = 

FTAN 

AMY(N,J)  ■ 

BMY(N,J)  i 
330   CONTINUE 

CMY(N,MC)  i 


Y  -  VELOCITY  INDUCED  AT  N-TH  CORE  VORTEX  BY  UNIT 
STRENGTH  DISTRIBUTED  SOURCE  ON  THE  J-TH  PANEL 

Y  -  VELOCITY  INDUCED  AT  N-TH  CORE  VORTEX  BY  UNIT 
STRENGTH  DISTRIBUTED  VORTEX  ON  THE  J-TH  PANEL 

.GT.  0)  GO  TO  530 
=  1,MM1 
XC(N" 
YC(N 

=  l,NODTOT 
XMID  -  X(J) 
XMID  -  X(J+1) 
YMID  -  Y(J) 
YMID  -  Y(J+1) 
5*ALOG( 


DX JP*DX JP+DY JP*DYJP ) / ( DXJ*DXJ+DYJ*DY J ) ) 
ATAN2(DYJP*DXJ-DXJP*DYJ,DXJP*DXJ+DYJP*DYJ) 
PI2INV*(FTAN*COSTHE(J)  -  FLOG*SINTHE ( J' 
PI2INV*(FLOG*COSTHE(J)  +  FTAN*SINTHE( J' 


=    1,MM1 
XC(N) 
YC(N) 
NP1 

XMID  - 
XMID  - 
YMID  - 
YMID  - 
\LOGt 


X{J) 

X(J+1) 
Y(J) 
"  >D 

:■:  JP*DX JP^-DY  JP*DY J?  )  /  ( DX J*DX J+DY J*DY J ) ) 


ATAN2 ( DYJP*DXJ-DXJP*DYJ , DXJP*DXJ+DY JP*DY J ) 
■  PI2INV*(FTAN*COSTHE(J)  -  FLOG*SINTHE( J) ) 
=  PI2INV*(FLOG*COSTHE(J)  +  FTAN*SINTHE ( J) ) 


CMX(N,MC) 


i  Y  -  VELOCITY  INDUCED  AT  N-TH  CORE  VORTEX  BY  UNIT 
STRENGTH  MC-TH  CORE  VORTEX  OTHER  THAN  ITSELF 

:  X  -  VELOCITY  INDUCED  AT  N-TH  CORE  VORTEX  BY  UNIT 
STRENGTH  MC-TH  CORE  VORTEX  OTHER  THAN  ITSELF 
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(NITR  .GT.  0)  RETURN 
N  =  1,MM1 
=  XC(N) 
=  YC(Nj 
MC  =  1,MM1 
.EQ.  N)  GO  TO  410 
XMID  -  XC(MC) 

:) 


=  YMID  -  YC(MC; 

=  SORT(DX*DX+DY*DY) 

=  DX/DIST 

=  DY/DIST 

=  -PI2INV*COSTHN/DIST 
=  +PI2INV*SINTHN/DIST 


IF 

DO  420 

XMID 

YMID 

DO  410 

IF  (MC 

DX 

DY 

DIST 

COSTHN 

SINTHN  = 

CMY(N,MC) 

CMX(N,MC) 
410   CONTINUE 
420   CONTINUE 

RETURN 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C     SUBROUTINE  KUTTA  (ALPHA, SINALF , COSALF, OMEGA, UX,UY)  C 

C  C 

C  USING  KUTTA  CONDITION  TO  DETERMINE  VORTICITY  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  KUTTA  (ALPHA, SINALF , COSALF , OMEGA, UX ,UY) 
COMMON  /BOD/  IFLAG,NLOWER,NUP?ER,NODTOT ,X(202 ) , Y(202) , 
+  COSTHE(201) ,SINTHE(201)/SS,NP1,MP2 

/COF/  A(201,211),NEQS 
/SING/  Q ( 200 ), GAMMA, QK( 200) ,GAMK 
/WAK/  VYW,VXW.WAKE,DT 
/CORV/  CV(20Q) ,XC(200),YC(200) ,M,TD , CCVX(200 ) ,CCVY(200) 


c 
c 
c 


c 
c 
c 


50 


120 


C 

c 
c 


c 
c 


110 
100 
130 


COMMON 
COMMON 
COMMON 
COMMON 
COMMON 
COMMON 
COMMON 


( 


/INF1/  AAN(201,201) ,BBN(201 ,201) ,AYNP1(201 
/INF2/  CCN(201,200) , CCT(201 ,200) ,CYNP1(200) ,CXN?i 
/GUST/  UG(200) ,VG(200) ,XGF,UGUST, VGUST 


3YNP1(201) 
(200) 


DIMENSION  B1(200),B2(200),AA(2),BB(2) 


DO  50 

B1(I) 
B2  (I) 


DO  130 
IF  (K 
IF  (K 
XMID 
YMID 
VTANG 


RETRIEVE  SOLUTION  FROM  A-MATRIX 

I  =  l,NODTOT 
=  A(I,NP1) 
=  A(I,NP2) 

FIND  VT  AT  TRAILING  EDGE  PANELS 


K  =  1 

.EQ.  1 
.EQ.  2 
=  0.5 
=  0.5 


,2 

I 

I 


=  1 

=  NODTOT 


( ( 1 . +UG ( I ) ) *COSALF-VG ( I ) *SINALF+UX ) *COSTHE ( I 
+  +  ((l.+UG(I))*SINALF+VG(I)*COSALF+UY)*SINTHE(I 

+  +  OMEGA*(YMID*COSTHE(I)  -  XMID^SINTHE (I) ) 

AA(K)  =  -  AAN(I,NP1)*SS/WAKE 

BB(K)  =  VTANG  +  AAN(I  ,NPl)*SS'lr GAMMA/ WAKE 

DO  120  J  =  1, NODTOT 

AA(K)  =  AA(K)  +  AAN(I,J)  -  BBN(I , J)*B1( J) 

3B(K)  =  3B(K)  -  3BN(I.J)*B2(J) 
CONTINUE 

ADD  CORE  VORTEX  CONTRIBUTION 


IF  (M 

MM1 

DO  110 

BB(K) 
CONTINUE 
CONTINUE 
CONTINUE 


.EQ.  1)  GO  TO  100 
=  M  -  1 
N  =  1,MM1 

BB(K)  +  CCI(I,N)*CW(U) 


SATISFYING  KUTTA  CONDITION 


SOLVE  FOR  VORTEX  STRENGTH 
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EE      =  AA(1)*AA(1)  -  AA(2)*AA(2) 

FF      =  AA(l)*BB(l)  -  AA(2)*BB(2)  -  SS/DT 

GG      =  BB(1)*BB(1)  "  BB(2)*BB(2)  +  2 . *SS*GAMMA/DT 

RADI     =  SQRT(FF*FF-EE*GG) 

GAMK    =  (-FF  -  RADI)/EE 
C 

C  CALCULATE  SOURCE  STRENGTH 

C 

DO  160   I  =  l.NODTOT 
160   QK(I)    =  GAMK^Bl(I)  +  B2(I) 

RETURN 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C      SUBROUTINE  NACA45 (Z , THICK, CAMBER, 3ETA)  C 

C  C 

C  EVALUATE  THICKNESS  AND  CAMBER  C 

C  FOR  NACA  4-  OR  5-DIGIT  AIRFOIL  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  NACA45 (Z , THICK , CAMBER , BETA) 
COMMON  /PAR/  NACA,TAU,EPSMAX,PTMAX 
THICK   =0.0 

IF  (Z  .LT.  l.E-10)     GO  TO  100 

THICK   =  5.*TAU*(.2969*SORT(Z)  -  Z*(.i26  +  Z*(.3537 
+  -  Z*(.2843  -  Z*.1015)))) 

100  IF  (EPSMAX  .EQ.  0.0)  GO  TO  130 
IF  (NACA  .GT.  9999)  GO  TO  140 
IF  (Z  .GT.  PTMAX)  GO  TO  110 

CAMBER   =  EPSMAX/PTMAX/?TMAX*(2.*?TMAX  -  Z)*Z 
DCAMDX  =  2. *EPSMAX/PTMAX/ PTMAX* (PTMAX  -  Z) 
GO  TO  120 
110   CAMBER   =  EPSMAX/ (1 . -PTMAX)**2*(1 .  +  Z  -  2 . *?TMAX) *( 1 .  -  Z) 

DCAMDX  =  2. *EPSMAX/ (1. -PTMAX )**2*( PTMAX  -  Z) 
120   BETA    =  ATAN( DCAMDX) 

RETURN 
130   CAMBER   =0.0 
BETA    =0.0 
RETURN 
140   IF  (Z  .GT.  PTMAX)       GO  TO  150 
W       =  Z/ PTMAX 

CAMBER   =  EPSMAX*W*((W  -  3.)*W  +  3.  -  PTMAX) 
DCAMDX  =  EPSMAX*3.*W*(1.  -  W) /PTMAX 
GO  TO  120 
150   CAMBER   =  EPSMAX*(1.  -  Z) 
DCAMDX   =  -  EPSMAX 
GO  TO  120 
END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C     SUBROUTINE  PRESS  (SINALF, COS ALF , OMEGA, UX,UY)  C 

c  c 

C  COMPUTE  UNSTEADY  FLOW  PRESSURE  DISTRIBUTION  C 

C  AND  VELOCITY  POTENTIAL  AT  MID-POINTS  OF  PANELS       C 

c  : 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  PRESS  ( SINALF , COSALF , OMEGA . JX , JY) 

HON  /BOD/  IFLAG,NLOWER,NUPPER,NODTOT,X  202   Z  202), 

+  COSTHE(201) ,SINTHE(201) ,SS,N?1,NP2 

COMMON  /CPD/  CP(200) 

COMMON  /NUM/  PI,PI2INV 

COMMON  /SING/  0(200) , GAMMA, OK (200) , GAMK 

COMMON  /WAK/  VYW, VXW . WAKE , DT 

COMMON  /CORV/  CV(200) ,XC(200 ) , YC(200) ,M,TD , CCVX(200) , CCVY( 200 ) 

COMMON  /INF1/  AAN(201,201),BBN(201,201),AYNP1(201) ,BYNP1(201) 

COMMON  /INF2/  CCN(201 , 200 ) , CCT( 201 , 200 ) , CYNP1 (200 ) , CXNP1 (200 ) 

COMMON  /POT/  PHI (200) ,PHIK(200) 

COMMON  /GUST/  UG(200) , VG(200) , XGF , UGUST, VGUST 


COMMON  /EXTV/  UE(200. 


US 


c 

c 
c 


c 

c 
c 


120 


140 
150 


C 
C 
C 
C 
C 


130 


C 

c 
c 


40 
20 


WRITE  (6,1000) 
FIND  TANGENTIAL  VELOCITY  VT  AT  MID-POINT  OF  I-TH  PANEL 


DO  130 

XMID 

YMID 

DX 

DY 

DIST 

VSX 

VSY 

VS. 

VTANG 


VTFREE 
VTANG 
DO  120 
VTANG 
CONTINUE 


I  =  1,N0DT0T 


0.5 
0.5  ' 

(x(i- 

CY(I+1) 


■1) 


X(I 

Y(I 


X(I  +  1 

Y(I  +  1 


50RT(DX*DX+DY*DY) 
(l.+UG(I) )*COSALF-VG(I 
(l.+UG(I 


*SINALF 
*C0SALF 


+  OMEGA*YMID 
-  OMEGA*XMID 


*SINALF+VGj(I 

VSX*VSX  +  VSY*VSY 

( ( 1 . +UG { I)  ) *C0SALF-VG ( I ) *SINALF+UX ) *C0STHE (I) 
+  ((l.+UG(I))*SINALF+VG(I)*COSALF+UY)*SINTHE(I) 

+  OMEGA* (YMID*COSTHE(I)  -  XMID*SINTHE (I) ) 

VTANG 

VTANG  +  SS*(GAMMA-GAMK)*AAN(I,NP1)/WAKE 

=  1,N0DT0T 

VTANG  -  BBN(I,J)*0K(J)  +  AAN ( I , J ) *GAMK 


UX 
UY 


ADD  CORE  VORTEX  CONTRIBUTION 


IF  (M  .EO 
MM1      = 
DO  140   N 
VTANG   = 
CONTINUE 
CONTINUE 
PHIK(I)  = 
C?(I)    = 
UE  ( I ) 
CONTINUE 


,  1)  GO  TO  150 
K  -  1 

=  1,MM1 
VTANG  + 


CCT(I,N)*CV(N) 


(VTANG- VTFREE ) *DIST 
VS  -  VTANG*VTANG 
VTANG 


COMPUTE  DISTURBANCE  POTENTIAL  BY  LINE  INTEGRAL  OF  VELOCITY  FIELD 

INTEGRATION  FROM  UPSTREAM  (AT  INFINITY)  TO  THE  LEADING  EDGE 

MPHI  =  10  *  NLOWER 

PINK  =0.0 

XL  =0.0 

DO  30  L  =  1,NPHI 

FRACT  =  FLOAT(L)/FLOAT(NPHI) 

XLP  =  -10.0  *  (1.0  -  COS(0.5*PI*FRACT)) 

DELX  =  XL  -  XLP 

XMID  =  0.5*(XL+XLP)*COSALF 

YMID  =  0.5*(XL+XLP)*SINALF 

XL  =  XLP 

VELX  =  UGUST 


ADD  CONTRIBUTION  OF  J-TH  PANEL 


DO  20 
DXJ 

DXJP 

DYJ 

DYJP 

FLOG 

FTAN 

CALMTJ 

SALMTJ 

APY 

BPY 

IF  (J 

VELX 

GO  TO 
VELX 
CONTINUE 


=  1,NP1 
XMID  -  X 
XMID  -  X 
YMID  -  Y 
YMID  -  7 


J)   X 

J+l) 

J) 

J+l) 


EQ 


. 5*ALOG ( ( DX JP*DX JP+DY JP*DY JP ) / ( DX J*DXJ+DYJ*DY J ) ) 
ATAN2 (DYJP*DXJ-DXJP*DYJ , DXJP*DXJ+DYJP*DYJ) 

SINALF*SINTHE(J) 

COSALF*SINTHE(J) 

+  FLOG*SALMTJ) 

-  FTAN*SALMTJ) 


20 


■COSALF*COSTHE(J 
-SINALF*COSTHE(J)    + 
PI2INV*(FTAN*CALMTJ 
PI2INV*(FL0G*CALMTJ 
.    NP1)    GO   TO   40 
VELX      -    BPY*QK(J)    +GAMK*APY 


=  VELX      +   SS*APY*(GAMMA-GAMK)/WAKE 
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C        ADD  CORE  VORTEX  CONTRIBUTION 
C 

IF  (M  .EC.  1)  GO  TO  50 

MH1      =  M  -  I 

DO  60   N  =  :,mmi 

DX      =  XMID  -  XC(N) 
DY      =  YMID  -  YC(N) 
DIST    =  SORT(DX*DX+DY*DY) 
COSTHN   =  DX/DIST 
SINTHN  =  DY/DIST 

SALMTN   =  -SINALF*CCSTHN  +  COSALF*SINTHN 
CPT      =  -?I2INV*SALMTN/DIST 
60    VELX    =  VELX   +  C?T*CV(N) 
50    CONTINUE 

PINK    =  PINK  +  VELX  *  DELX 
30    CONTINUE 
C 

C      COMPUTE  DISTURBANCE  POTENTIAL  AT  MID-POINT  OF  I-TH  PANEL 
C 

C        LOWER  SURFACE 
C 

DO  230   I  =  l,NLOWER 
PH      =  -PINK 
DO  240   J  =  I,MLCWER 
240   PH      =  PH  -  ?HIK(J) 

PHIK(I)  =  PH 
230   CONTINUE 

DO  270   I  =  l.NLOWER-1 
PHIK(I)  =  0.5*(PHIK(I)  +  PHIK(I+1)) 
270   CONTINUE 

PHIK(NLOWER)  =  0 . 5*(?HIK(NLOWER)  +  PINK) 
C 

C        UPPER  SURFACE 
C 

DO  250   I  =  MODTOT,NLCWER+1,-1 
PH      =  -PINK 
DO  260   J  =  NLOWER+l.I 
260   PH      =  PH  +  PHIK(J) 

?HIK(I)  =  PH 
250   CONTINUE 

DO  280   I  =  NODTOT,NLOWER+2,-l 
PHIK(I)  =  0.5*(PHIK(I)  +  PHIK(I-l)) 
280   CONTINUE 

PHIK(NLOWER+l)  =  0 .  5*(PHIK(NLOWER+l )  +  PINK) 
C 

C        COMPUTE  CP  AT  MID  POINT  OF  I-TH  PANEL 
C 

DO  290   I  =  l,NODTOT 
XMID     =  .5*(X(I)  +  X(I+1)) 
CP(I)    =  CP(I)  -  2.*(PHIK(I)-PHI(I))/DT 
WRITE  (6,1050)  I,XMID/QK(I)/GAMK,CP(I)/UE(I) 
290   CONTINUE 
1000  FORMAT (/,4X, ' J ' ,4X, ' X( J) ' , 6X, ' Q( J) ' , 5X, 'GAMMA' ,5X, 

+  'CP(J)1 /6X,'V(J)1 ,/) 
1050  FORMAT(I5,6F10.6) 
RETURN 
END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  SUBROUTINE  SETUP  "  C 

C  C 

C  SETUP  COORDINATES  OF  PANEL  NODES  AND  SLOPES  OF  PANELS  C 

C  COORDINATES  ARE  READ  FROM  INPUT  DATA  FILE  UNLESS  C 

C  THE  AIRFOIL  IS  OF  NACA  XXXX  OR  NACA  230XX  TYPE  C 

C  C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  SETUP 

COMMON  /BOD/  IFLAG, NLOWER ,MUPPER,NODTOT , X(202) , Y(202) , 
+        COSTHE(201) ,SINTHE(201) ,SS,NP1/NP2 
COMMON  /NUM/  PI,PI2INV 
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PI      =  3.1415926585 

c 
c 
c 

PI2INV  =  .5/PI 

SET  COORDINATES  OF  NODES  ON 

IF  (IFLAG  .NE.  0)  GO  TO  10 

NPOINT   =  NLOWER 

SIGN    =  -1.0 

NSTART   =  0 

DO  110   NSURF  =1,2 

DO  100  N  =  1, NPOINT 

FRACT    =  FLOAT(N-l)/FLOAT(NPOINT) 

Z       =  .5*(1.  -  COS(PI*FRACT)) 

I       =  NSTART  +  N 

CALL  BODY(Z,SIGN,X(I),Y(I)) 

100 

CONTINUE 

NPOINT   =  NUPPER 

SIGN    =1.0 

NSTART   =  NLOWER 

110 

CONTINUE 

NODTOT   =  NLOWER  +  NUPPER 

BODY  SURFACE 


X(NODTOT+l)  =  X(l) 
Y(NODTOT+l)  =  Y(l) 
GO  TO  20 
10    NODTOT   =  NLOWER  +  NUPPER 

READ  (1,501)  (X(I) .I=l,NODTOT+l) 
WRITE  (6,501)  (X(I),I=l,NODTOT+l) 
READ  (1,501)  (Y(I).I=l,NODTOT+l) 
WRITE  (6,501)  (Y(I),I=l,NODTOT+l) 


501 
20 

C 

C 

c 

FORMAT 
NPl 
NP2 

(6F10.6) 

=  NODTOT  +  1 
=  NODTOT  +  2 

SET  SLOPES  OF  PANELS 

SS 

DO  200 

DX 

DY 

DIST 

SS 

=  0.0 

I  =  1, NODTOT 
=  X(I+1)  -  xm 

=  y(i+i)  -  y(i)  - 
=  sqrt(dx*dx  +dy*dy) 

=  SS  +  DIST 

AND  CALCULATE  AIRFOIL  PERIMETER 


SINTHE(I)  =  DY/DIST 
COSTHE(I)  =  DX/DIST 
200   CONTINUE 
RETURN 
END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C     SUBROUTINE  TEWAK  (SINALF, COSALF)  C 

c  c 

C  COMPUTE  WAKE  ELEMENT  AT  THE  TRAILING  EDGE  C 

C  C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  TEWAK  (SINALF , COSALF) 

COMMON  /BOD/  IFLAG, NLOWER, NUPPER, NODTOT ,X(202) ,Y(202) , 
+  COSTHE(201) ,SINTHE(201) ,SS,NP1,NP2 

COMMON  /COF/  A( 20 1,211) ,NEOS 
COMMON  / S ING/  0(200), GAMMA  TOK ( 2 0 0 ) , GAMK 
COMMON  /WAK/  VYW , 7XW , WAKE , DT 
COMMON  /WAK2/  VYWK,VXWK 

COMMON  /CORV/  CV(200) ,XC(200) , YC(200 ) ,M,TD ,CCVX(200) , CCVY(200) 
COMMON  /INF1/  AAN(201/201),BBN(201/201),AYNP1(201),BYNP1(201) 
COMMON  /INF2/  CCN(201 , 200) , CCT(201 , 200 ) , CYNP1 (200) , CXNP1 (200 ) 
COMMON  /GUST/  UG(200) , VG(200) ,XGF ,UGUST, VGUST 
XMID    =  0.5  *  (X(NP1)  +  X(NP2)' 
YMID    =  0.5  *    (Y(NP1)  +  Y(NP2) 
UGW     =0.0 
VGW     =0.0 

XG      =  XMID*COSALF  +  YMID*SINALF 
IF  (XG  .GT.  XGF)  GO  TO  10 
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10 


c 

c 
c 
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UGW     =  UGUST 

VGW     =  VGUST 
VYWK    =  (l.+UGW)*SINALF+VGW*COSALF 

VXWK    =  (l.+UGW)*COSALF-VGW*SINALF 

DO  120   J  =  l,NODTOT 

VYWK    =  VYWK  +  AYN?1(J)*QK(J)  +  BYNP1 ( J)*GAMK 
VXWK    =  VXWK  -  BYNP1(J)*QK(J)  +  AYNP1 ( J)*GAMK 

ADD  CORE  VORTEX  CONTRIBUTION 


IF  (M  .EQ.  1)  GO  TO  140 

MM1      =  M  -  1 

DO  130  N  =  1,MM1 

VYWK    =  VYWK  +  CYNP1(N)*CV(N) 
130   VXWK    =  VXWK  +  CXNP1(N)*CV(N) 
140   CONTINUE 

RETURN 

END 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

C     SUBROUTINE  VELDIS (SINALF , COSALF)  C 

C  C 

C  COMPUTE  STEADY  FLOW  PRESSURE  DISTRIBUTION  C 

C  AND  VELOCITY  POTENTIAL  AT  MID-POINTS  OF  PANELS      C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCZCCCCCCCCCCCCC 

SUBROUTINE  VELDIS (SINALF, COSALF) 

COMMON  /BOD/  IFLAG ,NLOWER,NUPPER,NODTOT , X(202) , Y(202) , 
+  COSTHE(201) ,SINTHE(201) ,SS,NP1,NP2 

COMMON  /COF/  A (20 1,211) ,KUTTA 

COMMON  /CRD/  C?(2C0) 

COMMON  /NUM/  PI,PI2INV 

COMMON  /SING/  0(200) .GAMMA , OK (200) ,GAMK 

COMMON  /POT/  PHI(200),PHIK(2"00) 

COMMON  /GUST/  UG ( 200 ) ,VG( 200 ) ,XGF , UGUST, VGUST 

COMMON  /EXTV/  UE(200) 

WRITE  (6,1000) 


C 
C 
C 


50 


C 

c 

c 


c 
c 
c 


RETRIEVE  SOLUTION  FROM  A-MATRIX 

DO  50    I  =  l,NODTOT 
0(1)     =  A(I,KUTTA+1) 
GAMMA    =  A(KUTTA,KUTTA+1) 

FIND  VT  AND  CP  AT  MID-POINT  OF  I-TH  PANEL 

DO  130  I  =  l,NODTOT 

XMID  =  .5*(X(I)  +  X(I+1)) 

YMID  =  .5*(Y(I)  +  Y(I+1)) 

VTANG  =  COSALF*COSTHE(I)  +  SINALF*SINTHE ( I ) 

VTFREE  =  VTANG 

ADD  CONTRIBUTION  OF  J-TH  PANEL 

DO  120   J  =  l,NODTOT 

FLOG    =0.0 

77  AN     =  PI 

IF  ;j  .EQ.  I   30  TO  LOO 


100 


120 


LJXJ 

DXJP 

DYJ 

DYJP 

FLOG 

FTAN 
CTIMTJ 

STIMTJ 

AA 

B 

VTANG 
CONTINUE 


XMID  - 
XMID  - 
YMID  - 
YMID  - 


A(J) 


X(J+1) 

Y  J) 

Y(J+1) 
5*ALOG( (DXJP*DXJP+DYJP*DYJP ) / (DXJ*DXJ+DYJ*DYJ ) ) 
=  ATAN2(DYJP*DXJ-DXJP*DYJ,DXJP*DXJ+DYJP*DYJ) 
=  COSTHE(I)*COSTHE(J)  +  SINTHE ( I )*SINTHE ( J) 
=  SINTHE(I)*COSTHE(J)  -  COSTHE (I ) ^SINTHE ( J) 
=  PI2INV*(FTAN'lrCTIMTJ  +  FLOG*STIMTJ) 
=  PI2INV*(FLOG*CTIMTJ  -  FTAN*STIMTJ) 
=  VTANG  -  B*Q(J)  +GAMMA*AA 
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CP(I)    =1.0-  VTANG*VTANG 

US(I)    =  VTANG 

WRITE  (6,1050)  I,XMID,Q(I),GAMMA,CP(I),UE(I) 
C 

C        INITIAL  SET-UP  "OR  DISTURBANCE  POTENTIAL  CALCULATION 
C 

DX      =  X(I+1)  -  X(I) 

DY      =  7(1+1)  -  Y(I) 

DIST    =  SQRT'DX*DX+DY*DY) 

PHI  (I)   =  (VTANG-VTFREE)*DIST 
130   CONTINUE 
C 

C     COMPUTE  DISTURBANCE  POTENTIAL  BY  LINE  INTEGRAL  OF  VELOCITY  FIELD 
C 

C        INTEGRATION  FROM  UPSTREAM  (AT  INFINITY)  TO  THE  LEADING  EDGE 
C 

NPHI     =  10  *  NLOWER 

PIN     =0.0 

XL      =0.0 

DO  30    L  =  1,MPHI 

FRACT    =  FLOAT (L) /FLOAT (NPHI) 

XLP     =  -10.0  *  (1.0  -  COS(0.5*PI*FRACT)) 

DELX    =  XL  -  XLP 

XMID     =  0.5*(XL+XL?)*COSALF 

YMID     =  0.5*(XL+XLP)*SINALF 

XL      =  XLP 

VELX    =   UGUST 
C 

C  ADD  CONTRIBUTION  OF  J-TH  PANEL 

C 

DO  20    J  =  l,NODTOT 

DXJ     =  XMID  -  X(J) 

DXJP     =  XMID  -  X(J+1) 

DYJ     =  YMID  -  Y(J) 

DYJ?    =  YMID  -  Y(J+1) 

FLOG    =  .5*ALOG((DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ)) 

FTAN    =  ATAN2(DYJP*DXJ-DXJP*DYJ,DXJP*DXJ+DYJP*DYJ) 

CALMTJ   =  -COSALF*COSTHE(J)  -  SINALF*SINTHE ( J) 

SALMTJ   =  -SINALF*COSTHE(J)  +  COSALF*SINTHE ( J) 

APY     =  ?I2INV*(FTAN*CALMTJ  +  FLOG*SALMTJ) 

BPY     =  ?I2INV*(FLOG*CALMTJ  -  FTAN*SALMTJ) 

VELX    =  VELX   -  BPY*Q(J)  +GAMMA*APY 
20   CONTINUE 

PIN     =  PIN  +  VELX  *  DELX 
30    CONTINUE 
C 

C      COMPUTE  DISTURBANCE  POTENTIAL  AT  MID-POINT  OF  I-TH  PANEL 
C 

C        LOWER  SURFACE 
C 

DO  230   I  =  1, NLOWER 

PH      =  -PIN 

DO  240   J  =  I, NLOWER 
240   PH      =  PH  -  PHI (J) 

PHI (I)   =  PH 
230   CONTINUE 

DO  270   I  =  1. NLOWER- 1 

PHI(I)   =  0.5*' PHI (I)  +  ?HI(I+1)) 
270   CONTINUE 

PHI (NLOWER)  =  0.5* (PHI (NLOWER)  +  PIN) 
C 

C        UPPER  SURFACE 
C 

DO  250   I  =  NODTOT,NLOWER+l,-l 

PH      =  -PIN 

DO  260   J  =  NLOWER+1,1 
260   PH      =  PH  +  PHI (J) 

PHI (I)   =  PH 
250   CONTINUE 

DO  280   I  =  NODTOT,NLOWER+2,-l 
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PHI(I)   =  0.5*(?HI(I)  +  PHI(I-l)) 

230   CONTINUE 

PHI'NLOWER+1)  =  0 . 5*(PHI (NLOWER+1 )  +  PIN) 
1000  F0RMAT(/,4X,  '  J1  ,4X,  'X(J)  '  ,6X,  'Q(J)  '  ,5X,  'GAMMA'  ,5X, 

+   C?IJ)',5X,'V(J)',/) 
1050  ?ORMAT(I5,5FI0.6) 

RETURN 

END 
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APPENDIX  B 
EXAMPLE  INPUT  DATA  FOR  PROGRAM  U2DIIF 


ll 


THIS  IS  AN  EXAMPLE  OUTPUT  DATA  OBTAINABLE  FROM  PROGRAM  U2DIIF 

AIRFOIL  :  MI5ES  8.4%  THICKNESS  (COORDINATES  ARE  INPUT  BY  USER) 

PANEL  NUMBER  :  NLOWER  =  25  ,  NUPPER  =25 

AIRFOIL  MOTION  :  MODIFIED  RAMP  AOA  CHANGE  ABOUT  MID  CHORD 

INITIAL  AOA  :  2.5  DEGREES 

FINAL  AOA  :  7.5  DEGREES 

AOA  RISE  TIME  :  1.5  CHORD  LENGTHS 

COMPUTATION  TIME  STEP  :  0.05  DURING  TRANSIENT  MOTION,  INCREASES 

PROGRESSIVELY  AFTER  FINAL  AOA  IS  REACHED 
*********2*********2*********3*********4**^**^^**5*^*^^^*^*6^*x^*^^;ic'i:7 


01   25 

1.000000 

0 

25 

.994858 

0 

.980866 

0 

.958884 

0 

.929536 

0 

.893455 

0.851308 

0 

.803815 

0 

.751753 

0 

.695943 

0 

.637271 

0 

.576620 

0.514918 

0 

.453098 

0 

.392084 

0 

.332794 

0 

.276105 

0 

.222365 

0.173361 

0 

.129819 

0 

.091393 

0 

.059146 

0 

.033560 

0 

.015010 

0.003767 

0 

.000000 

0 

.003767 

0 

.015003 

0 

.033560 

0 

.059146 

0.091393 

0 

.129819 

0 

.173861 

0 

.222865 

0 

.276105 

0 

.332791 

0.392082 

0 

.453095 

0 

.514915 

0 

.576617 

0 

.637266 

0 

.695946 

0.751750 

0 

.803815 

0 

.851308 

0 

.893455 

0 

.929536 

0 

.958884 

0.980866 

0 

.994858 

1 

.000000 

0.000000 

-o. 

.000732 

-0 

.002784 

-0 

.005721 

-0 

.009351 

-0 

.013459 

-0.017837 

-0. 

022235 

-0. 

026618 

-0. 

030671 

-0. 

034289 

-0. 

037341 

-0.039712 

-0. 

041314 

-0. 

042083 

-0. 

041979 

-0. 

,040979 

-0. 

039096 

-0.036360 

-0. 

032820 

-0. 

028555 

-0. 

023651 

-0. 

018220 

-0. 

012379 

-0.006259 

0. 

000000 

0. 

006259 

0. 

012379 

0. 

018220 

0. 

023651 

0.028555 

0 

.032820 

0 

.036360 

0 

.039096 

0 

.040979 

0 

.041979 

0.042083 

0 

.041314 

0 

.039712 

0 

.037341 

0 

.034289 

0 

.030671 

0.026618 

0 

.022285 

0 

.017837 

0 

.013459 

0 

.009351 

0 

.005721 

0.002784 

0 

.000782 

0 

.000000 

2.50000 

5 

.000000 

1.5 

0.0 

0.5 

0.0 

0.000000 

0 

.000000 

0 

.000000 

2.000000 

0 

.050000 

0.0001 

0 

.000000 

0.00 
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APPENDIX  C 
EXAMPLE  OUTPUT  DATA  FROM  PROGRAM  U2DIIF 


DATA  READ  FROM  FILE  CODE  1 

11 

THIS  IS  AN  EXAMPLE  OUTPUT  DATA  OBTAINABLE  FROM  PROGRAM  U2DIIF 

AIRFOIL  :  MISES  8.4%  THICKNESS  (COORDINATES  ARE  INPUT  BY  USER) 

PANEL  NUMBER  :  NLOWER  =  25  ,  NUPPER  =  25 

AIRFOIL  MOTION  :  MODIFIED  RAMP  AOA  CHANGE  ABOUT  MID  CHORD 

INITIAL  AOA  :  2.5  DEGREES 

FINAL  AOA  :  7.5  DEGREES 

AOA  RISE  TIME  :  1.5  CHORD  LENGTHS 

COMPUTATION  TIME  STEP  :  0.05  DURING  TRANSIENT  MOTION,  INCREASES 

PROGRESSIVELY  AFTER  FINAL  AOA  IS  REACHED 
^^^A^A^^yc]_^^A*^^^*^2^^^*^^^*^3**^^^^^^^4^*^^^^^^^5*^x^^^^^^6^^^^^^'^,ir'^7 


1   25 
1.000000 

25 
0.994858 

0.980866 

0 

.958884 

0 

.929536 

0 

.393455 

0.851308 

0.803815 

0.751753 

0 

.695948 

0 

.637271 

0 

.576620 

0.514918 

0.453098 

0.392084 

0 

.332794 

0 

.276105 

0 

.222855 

0.173361 

0.129819 

0.091393 

0 

.059146 

0 

.033560 

0 

.015010 

0.003767 

0.000000 

0.003767 

0 

.015008 

0 

.033560 

0 

.059146 

0.091393 

0.129319 

0.173861 

0 

.222865 

0 

.275105 

0 

.332791 

0.392082 

0.453095 

0.514915 

0 

.576617 

0 

.637266 

0 

.695946 

0.751750 

0.803815 

0.851308 

0 

.893455 

0 

.929536 

0 

.958884 

0.980866 

0.994858 

1.000000 

0.000000 

-0.000782 

-0.002784 

-0 

.005721 

-0 

.009351 

-0 

.013459 

0.017837 

-0.022285 

-0.026618 

-0. 

030671 

-0. 

034289 

-0. 

037341 

0.039712 

-0.041314 

-0.042083 

-0. 

041979 

-0. 

040979 

-0. 

039096 

0.036360 

-0.032820 

-0.028555 

-0. 

023651 

-0. 

018220 

-0. 

012379 

0.006259 

0.000000 

0.006259 

0. 

012379 

0. 

018220 

0. 

023651 

0.028555 

0.032820 

0.036360 

0 

.039096 

0 

.040979 

0 

.041979 

0.042083 

0.041314 

0.039712 

0 

.037341 

0 

.034289 

0 

.030671 

0.026618 

0.022285 

0.017837 

0 

.013459 

0 

.009351 

0 

.005721 

0.002784 

0.000782 

0.000000 

2.500000 

5.000000 

1.500000 

0 

.000000 

0 

.500000 

0 

.000000 

0.000000 

0.000000 

0.000000 

0.000000 

2.000000 

0.050000 

0.000100 

0 

.000000 

AIRFOIL  PERIMETER  LENGTH  ■   2.018599 


STEADY  FLOW  SOLUTION  AT  ALPHA  =    2.500000 


X(J) 


Q(J) 


GAMMA 


:f<j) 


7(J) 


0 

.997429 

0 

.355723 

j 

.074003 

0 

.316305 

-0 

.325359 

2 

0 

.987862 

0 

.356105 

0 

.074003 

0 

.206074 

-0 

.691025 

3 

0 

.969875 

0 

.365026 

0 

.074003 

0 

.133790 

-0 

.930704 

4 

0 

.944210 

0 

.378836 

0 

.074003 

0 

.082276 

-0 

.957979 

5 

0 

.911495 

0 

.394973 

0 

.074003 

0 

.043033 

-0 

.978247 

6 

0 

.872381 

0 

.412926 

0 

.074003 

0 

.012034 

-0 

.993965 

7 

0 

.827561 

0 

.432568 

0 

.074003 

-0 

.012724 

-1 

.006342 

8 

0 

.777784 

0 

.453710 

0 

.074003 

-0 

.032414 

-1 

.016078 

9 

0 

.723850 

0 

.476112 

0 

.074003 

-0 

.048010 

-1 

.023724 

10 

0. 

666609 

0. 

500047 

0. 

074003 

-0. 

059905 

-1. 

029517 

11 

0. 

606945 

0. 

525455 

0. 

074003 

-0. 

068405 

-1. 

033637 

12 

0. 

,545769 

0. 

552654 

0. 

074003 

-0. 

073563 

-1. 

036129 
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13 

0, 

.484008 

0, 

.531715 

0 

.074003 

-0, 

.075309 

-1, 

.036971 

14 

0, 

.422591 

0, 

.612382 

0 

.074003 

-0. 

.073531 

-]_  ( 

.036114 

15 

0. 

.362439 

0, 

.646430 

0, 

.074003 

-0. 

,067923 

-l! 

.033406 

16 

0, 

,304449 

0, 

.633297 

0 

.074003 

-0, 

,057653 

-l 

.023430 

17 

0. 

.249435 

0, 

.723595 

0 

.074003 

-0, 

.041391 

-i 

.020731 

18 

0, 

,198363 

0. 

.  763523 

0, 

.074003 

-0, 

,019114 

-i, 

.009512 

19 

0, 

.151340 

0. 

.319732 

0 

.074003 

0, 

,013374 

-0, 

.993290 

20 

0, 

.110606 

0, 

.379132 

0 

.074003 

0. 

,059669 

-0, 

.969707 

21 

0 

.075269 

0. 

.951277 

0 

.074003 

0. 

.127774 

-0, 

.933931 

22 

0, 

.046353 

1 
J-  1 

.043426 

0 

..074003 

0, 

.232984 

-0, 

.375795 

23 

0 

.024285 

1 

1 

.172734 

0 

.074003 

0 

.409236 

-0 

.763612 

24 

0, 

.009383 

1. 

.380641 

0 

.074003 

0 

.727179 

-0 

.522323 

25 

0, 

.001384 

1. 

,653644 

0 

.074003 

0, 

,945112 

0 

.234282 

26 

0, 

,001884 

0, 

.545367 

0 

.074003 

-0, 

.815326 

1 

.347341 

27 

0. 

.009387 

-0. 

,275210 

0 

.074003 

-1, 

.084184 

1 

.443670 

23 

0, 

.024234 

-0, 

,497235 

0 

.074003 

-0 

.372601 

1 

.363430 

29 

0, 

.046353 

-0, 

,580394 

0 

.074003 

-0 

.723499 

1 

.312821 

30 

0, 

.075269 

-0, 

.613032 

0 

.074003 

-0, 

.624167 

1 

.274428 

31 

0, 

.110606 

-0, 

.526699 

0 

.074003 

-0 

.552954 

1 

.246176 

32 

0, 

.151840 

-0, 

.645594 

0 

.074003 

-0, 

.498371 

1 

.224080 

33 

0, 

,193363 

-0, 

.649755 

0 

.074003 

-0. 

.453794 

1 

.205734 

34 

0, 

,249485 

-0, 

.650971 

0 

.074003 

-0, 

.415592 

1 

.189787 

35 

0, 

,304448 

-0. 

.650596 

0 

.074003 

-0. 

.381557 

1 

.175397 

36 

0 

,362436 

-0 

.649624 

0 

.074003 

-0 

.349957 

1 

.161877 

37 

0, 

.422538 

-0. 

.648020 

0 

.074003 

-0, 

.319875 

1 

.148858 

38 

0, 

.434005 

-0 

.646281 

0 

.074003 

-0 

.290579 

1 

.136037 

39 

0, 

,545766 

-0 

.544572 

0 

.074003 

-0 

.261352 

1 

.123099 

40 

0 

,606941 

-0. 

.642990 

0 

.074003 

-0 

.231604 

1 

.109776 

41 

0, 

.  666606 

-0, 

.541396 

0 

.074003 

-0, 

.200941 

1 

.095875 

42 

0 

.723848 

-0 

.639980 

0 

.074003 

-0, 

.168763 

1 

.081094 

43 

0, 

.777782 

-0, 

.633504 

0 

.074003 

-0 

.134640 

1 

.065195 

44 

0 

.327561 

-0, 

.537214 

0 

.074003 

-0 

.097676 

1 

.047701 

45 

0 

.372381 

-0, 

.535380 

0 

.074003 

-0 

.056864 

1 

.023039 

46 

0 

.911495 

-0 

.634260 

0 

.074003 

-0 

.010900 

1 

.005436 

47 

0 

.944210 

-0 

.632196 

0 

.074003 

0 

.042413 

0 

.978564 

48 

0 

.969875 

-0 

.629793 

0 

.074003 

0 

.107609 

0 

.944665 

49 

0 

.987362 

-0. 

.625960 

0 

.074003 

0 

.193062 

0 

.898297 

50 

0 

.997429 

-0 

.619334 

0 

.074003 

0 

.316307 

0 

.826858 

CD 

= 

0.000829 

CL  = 

•     1 

3.303076 

CM  =  - 

■0.' 

080325 

*************************************** 

***   BEGIN  UNSTEADY  FLOW  SOLUTION   **** 
*************************************** 


TIME  STEP  TK  =   0.050000 


TK  -  TKM1  = 


0.050000 


ALPHA (T)  = 
U(T) 


2.516295 
0.000000 


OMEGA(T)  = 
V(T) 


-0.011248 

-0.005624 


NITR 

0 
1 
2 
3 


VXW 

0.999048 
0.907832 
0.904138 
0.903985 


7YW 

0.043619 
0.005991 
0.007297 
0.007241 


WAKE 

0.050000 
0.045393 
0.045208 
0.045201 


THETA 

0.043633 
0.006600 
0.008070 
0.008010 


GAMMA 

0.740032E-01 
0.744799E-01 
0.744662E-01 
0.744652E-01 


CONVERGED  SOLUTION  OBTAINED  AFTER  NITR  =    3 

J    X(J)       Q(J)      GAMMA     CP(J)       V(J) 
1   0.997429   0.435837   0.074466   0.299470  -0.838202 
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1 

4 

5 
6 

7 

3 

9 

10 

:: 

12 
13 
14 
15 

16 

i  -t 

13 
19 
20 
21 
22 
23 
24 
25 
26 
27 
28 
29 
30 
31 
32 
33 
34 
35 
36 
37 
33 
39 
40 
41 
42 
43 
44 
45 
46 
47 
43 
49 
50 


.987362 
.969375 
.944210 
.911495 
.372381 
.327561 
.777784 
.722850 
.666609 
.606945 
.545769 
.484003 
.422591 
.362439 
.304449 
.2^9485 
.193363 
.151340 
.110606 
.075269 
.046353 
.024235 
.009338 
.001384 
.001334 
.009387 
.324234 
.046353 
.075269 
.110606 
.151340 
.198263 
.249435 
.304443 
.362426 
.422588 
.484005 
.545766 
.606941 
.  666606 
.722843 
.777782 
.827561 
.872381 
.911495 
.944210 
.969375 
.987862 
.997429 


0.422S52 

0.421793 
0.427256 
0.435974 
0.447131 
0.460598 
0.475913 
0.492820 
0.511571 
0.532044 
0.554564 
0.579201 
606206 
635915 
669133 
706140 
748100 
796633 
353325 
924099 
014820 
143358 
351349 
534712 
564272 
246433 
467313 
551795 
590860 
611392 
622549 
529333 
533515 
536433 
539059 
0.641343 
0.643768 
0.646433 
649578 
652913 
656699 
660707 
665236 
0.670135 
0.675261 
680614 
686543 
692719 
699982 


0 
0 
0 
0 
0 
0 
0 
0 

1 

1 


0 
-0 
-0 
-0 
-0 
-0 
-0 
-0 

-a 

-0 
-0 


-0 

•J 

-0 
-0 

-0 


•0 
•0 
•0 
■0 


0.074466 
0.074466 
0.074466 
0.074466 
0.074465 
0.074466 
.074466 
.074465 
074466 
074465 
074466 
,074466 
0.074466 
0.074466 
0.074465 
.074466 
074466 


0 

0 

0. 

0. 

0. 

0. 


0 
Q 
0 
0 
0 
0 


074465 
074466 
074466 
074466 
0.074466 
0.074465 
074466 
074465 
074466 
074466 
0.074466 
0.074466 
074466 
074466 
074466 
374466 
074466 
074466 
074466 
074466 
,074466 
0.074466 
0.074466 
0.074466 
0.074466 
0.074466 
0.074466 
0.074466 
0.074466 
0.074466 
0.074466 
0.074466 


0 
0 
0 

0 


0 
0 
0 
0 
0 
0 
0 
0 
0 


0.196120 

0.132077 

0.088439 

0.055859 

0.029322 

0.008149 

.010442 

.026773 

.041132 

.053494 

.063531 

.070990 

.075185 

.075487 

.070722 

.359745 

.040826 

.011154 

0.033396 

0.100715 

205876 

382654 

.703606 

952740 

746972 

037466 

338061 

■0.593563 

■0.596473 

■0.527117 

•0.474313 

■0.433321 

■0.399064 

■0.369826 

■0.343641 

■0.319346 

■0.295850 

■0.272083 

247068 


-0 
-0 
■0 
■0 
■0 
-0 
■0 
-0 
-0 
-0 
•0 
■0 


0 
0 
0 
0 

■0 

•1 

■0 


-1 
-1 
-1 


•0 
■0 
■0 
•0 
•0 


220048 


190134 
156552 
118313 
•0.074279 
•0.023233 
0.036802 
0.109850 
0.203356 
0.333160 


-0.899569 
-0.936976 
-0.962416 
-0.981131 
-0.995676 
-1.007037 
-1.015962 
-1.022942 
-1.028238 
032001 
034262 
035018 
-1.034204 
-1.031672 
-1.027008 
-1.019773 
-1.009171 
-0.993762 
-0.971231 
-0.936853 
-0.880676 
-0.776542 
535871 
209596 
322641 
.430099 
360476 
.307915 
271431 
244626 
223580 
206047 
190715 
176793 
163586 
150745 
137963 
124936 
111386 
.097122 
.081349 
065285 
046980 
.026307 
002476 
.974108 
938386 
.889792 
815602 


-0 
0 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
0 
0 
0 
0 


CD  =      0.001539 


CL  =      0.302054 


CM  =   -0.088450 


TRAILING  VORTICES    DATA 

M  X(M)  Y(M)  CIRC 

1       1.322599      3.300131    -0.000933 


TIME   STEP   TK  =        0.749999 


TK   -   TKM1   = 


0.050000 


ALPHA(T)    =        4.999996 
U(T)  =        0.000000 


OMEGA(T)    =      -0.087266 
V(T)  =      -0.043633 


NITR  VXW  VYW  WAKE  THETA  GAMMA 

0      0.905684   -0.000916      0.045284   -0.001012      0.103235E+00 
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1      0.905735   -0.000649      0.045287    -0.000717      0 .106563E+00 


CONVERGED   SOLUTION  OBTAINED   AFTER  NITR   = 


X(J) 


Q(J) 


GAMMA 


CP(J) 


v(J) 


1 

0 

.997429 

1 

.115997 

0 

.106565 

0 

.311649 

-0 

.908159 

2 

0 

.987862 

1 

.060333 

0 

.106565 

0 

.221106 

-0 

.957033 

3 

0 

.969875 

1 

.031653 

0 

.106565 

0 

.170306 

-0 

.983709 

4 

0 

.944210 

1 

1 

.013085 

0 

.106565 

0 

.141163 

-0 

.998976 

5 

0 

.911495 

0 

.998141 

0 

.106565 

0 

.123949 

-1 

.008098 

6 

0 

.872331 

0 

.985265 

0 

.106565 

0 

.114073 

-1 

.013416 

7 

0 

.827561 

0 

.974241 

0 

.106565 

0 

.109016 

-1 

.015142 

8 

0 

.777784 

0 

.964993 

0 

.106565 

0 

.107136 

-1 

.017013 

9 

0 

.723850 

0 

.957451 

0 

.106565 

0 

.107130 

-1 

.016594 

10 

0. 

,666609 

0, 

,952120 

0. 

,106565 

0, 

,108458 

-1, 

,015039 

11 

0, 

,506945 

0, 

,949235 

0, 

,106565 

0, 

,110657 

-1, 

,012669 

12 

0. 

,545769 

0, 

,949393 

0. 

,106565 

0, 

,113665 

-1, 

,009317 

13 

0. 

,434008 

0, 

,952960 

0. 

,106565 

0. 

,117540 

-1, 

,004971 

14 

0, 

,422591 

0, 

,960451 

0. 

,106565 

0, 

,122480 

-0, 

,999505 

15 

0. 

,362439 

0, 

,972455 

0. 

,106565 

0, 

,128967 

-0, 

,992654 

16 

0, 

,304449 

0, 

,990059 

0. 

,106565 

0. 

,138078 

-0, 

,983855 

17 

0, 

,249485 

1, 

,013807 

0. 

,106565 

0. 

,150962 

-0, 

,972486 

13 

0, 

,198363 

1, 

,045081 

0. 

,106565 

0, 

,169616 

-0, 

,957451 

19 

0, 

,151840 

1, 

,035785 

0. 

,106565 

0, 

,197364 

-0, 

,936848 

20 

0, 

,110606 

1, 

,133034 

0, 

,106565 

0, 

,239070 

-0, 

,907675 

21 

0, 

,075269 

1, 

,206453 

0. 

,106565 

0, 

,303822 

-0, 

,363894 

22 

0, 

,046353 

1, 

,298127 

0. 

,106565 

0, 

,407818 

-0, 

,793054 

23 

0, 

,024285 

1. 

,428849 

0, 

,106565 

0, 

,583190 

-0, 

,663132 

24 

0, 

,009388 

1, 

,631804 

0, 

,106565 

0, 

,872165 

-0, 

.369636 

25 

0, 

,001884 

1, 

,819807 

0, 

,106565 

0, 

,765241 

0, 

.485826 

26 

0, 

,001334 

0, 

,373179 

0, 

,106565 

-1, 

,559307 

1, 

.595823 

27 

0, 

,009387 

-0, 

.529374 

0, 

,106565 

-1, 

,564167 

1 

.590937 

28 

0, 

,024234 

-0, 

.755145 

0, 

,106565 

-1. 

,202155 

1 

.468063 

29 

0, 

,046353 

-0 

.836350 

0, 

,106565 

-0, 

,991312 

1 

.389589 

30 

0, 

,075269 

-0, 

.874103 

0, 

,106565 

-0, 

,864650 

1, 

.338450 

31 

0, 

.110606 

-0, 

.896239 

0, 

,106565 

-0, 

,781035 

1 

.302194 

32 

0, 

,151840 

-0, 

.912096 

0, 

,106565 

-0, 

.721021 

1, 

.274529 

33 

0, 

,198363 

-0, 

.926607 

0, 

,106565 

-0, 

,674005 

1, 

.251843 

34 

0, 

,249435 

-0, 

.941343 

0, 

,106565 

-0, 

,634257 

1, 

.232132 

35 

0, 

,304448 

-0, 

,957406 

0, 

,106565 

-0, 

,598321 

1, 

.214143 

36 

0, 

,362436 

-0, 

.975544 

0, 

,106565 

-0, 

,563539 

1, 

.196886 

37 

0, 

,422588 

-0, 

.995441 

0, 

,106565 

-0, 

,528592 

1, 

.179831 

38 

0, 

,484005 

-1, 

.017302 

0, 

,106565 

-0, 

,492391 

1 

.162518 

39 

0, 

,545766 

-1, 

.041010 

0, 

,106565 

-0, 

,454043 

1 

.144540 

40 

0, 

.606941 

-1, 

.066387 

0, 

.106565 

-0, 

,412924 

1, 

.125555 

41 

0, 

.666606 

-1, 

.093023 

0, 

.106565 

-0, 

.368740 

1, 

.105335 

42 

0, 

.723848 

-1 

.120818 

0, 

.106565 

-0, 

.321026 

1, 

.083543 

43 

0, 

.777782 

-1 

.149241 

0, 

.106565 

-0, 

.269645 

1 

.059939 

44 

0, 

.827561 

-1 

.178310 

0, 

.106565 

-0, 

.213999 

1 

.034027 

45 

0, 

.872381 

-1. 

.207644 

0, 

.106565 

-0, 

.153563 

1 

.005269 

46 

0, 

.911495 

-1 

.236895 

0, 

.106565 

-0, 

.087684 

0 

.972965 

47 

0, 

.944210 

-1 

.265976 

0, 

.106565 

-0, 

.014749 

0 

.935775 

48 

0, 

.969875 

-1, 

.296105 

0, 

.106565 

0, 

.069103 

0 

.890819 

49 

0, 

.987862 

-1 

.330154 

0, 

.106565 

0, 

.171240 

0 

.832327 

50 

0 

.997429 

-1 

.380203 

0 

.106565 

0 

.308161 

0 

.746084 

CD 

= 

0.033387 

r^r      — 
LL    - 

3.645338 

CM   =    ■ 

•0. 

224298 

TRAILING  VORTICES   DATA 
M  X(M)  Y(M) 


CIRC 


1 

1 

.700760 

0, 

.077052 

-0, 

.000933 

2 

1 

.651289 

0, 

.072156 

-0, 

.001625 

3 

1 

.602002 

0, 

.066331 

-0, 

.002229 

4 

1 

.552845 

0, 

.060031 

-0, 

.002784 

5 

1 

.503779 

0 

.053471 

-0 

.003304 

6 

1 

.454797 

0 

.046788 

-0 

.003797 

7 

1 

.405895 

0, 

.040107 

-0 

.004256 
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8  1.357057   0.033554  -0.004687 

9  1.308305   0.027216  -0.005100 


10 

1 

.259669 

0 

.021159 

-0, 

.005481 

11 

1 

.211178 

0 

.015541 

-0, 

.005801 

12 

1 

.162916 

0 

.010481 

-0, 

.006100 

13 

1 

.115048 

0, 

.006079 

-0, 

.006349 

14 

1 

.067918 

0, 

.002406 

-0, 

.006559 

15 

1 

.022643 

-0, 

.000016 

-0, 

.006720 

TIME  STEP  TK  =   1.449992 


TK  -  TKM1  =   0.050000 


ALPHA (T)  =   7.483698 
U(T)     =   0.000000 

OMEGA (T) 
V(T) 

=   -0.011249 
=   -0.005625 

NITR    VXW       VYW 

WAKE 

THETA     GAMMA 

0  0.901699   0.009872 

1  0.901200   0.011028 

0.045088 
0.045063 

0.010948   0.145377E+00 
0.012236   0.146997E+00 

CONVERGED  SOLUTION  OBTAINED  AFTER  NITR  = 


X(J) 


Q(J) 


GAMMA 


CP(J) 


v(J) 


1 

0 

.997429 

1 

.101898 

0 

.146996 

0 

.332408 

-0 

.364078 

2 

0 

.987862 

1 

.099609 

0 

.146996 

0 

.228004 

-0 

.921138 

3 

0 

.969875 

1 

.116333 

0 

.146996 

0 

.160643 

-0 

.955077 

A 

0 

.944210 

1 

.140985 

0 

.146996 

0 

.114863 

-0 

.976579 

i 

0 

.911495 

1 

.163619 

0 

.146996 

0 

.032627 

-0 

.990894 

6 

0 

.372381 

1 

.197904 

0 

.146996 

0 

.060523 

-1 

.000257 

7 

0 

.827561 

1 

.228671 

0 

.146996- 

0 

.046710 

-1 

.005378 

8 

0 

.777784 

1 

.260826 

0 

.146996 

0 

.039936 

_  1 

.008490 

9 

0 

.723850 

1 

.294203 

0 

.146996 

0 

.039112 

-1 

.008660 

10 

0, 

,666609 

1, 

,329222 

0, 

,146996 

0, 

,043766 

-1, 

,006560 

11 

0, 

,606945 

1, 

,366027 

0, 

,146996 

0, 

,053332 

-1. 

.002352 

12 

0. 

,545769 

1, 

,405159 

0, 

,146996 

0, 

,067653 

-0. 

,995947 

13 

0, 

,484008 

1, 

,446882 

0, 

,146996 

0, 

,086578 

-0, 

,987214 

14 

0. 

,422591 

1, 

,491619 

0, 

,146996 

0, 

,110103 

-0, 

.975912 

15 

0, 

,362439 

1, 

,539856 

0. 

,146996 

0, 

,138557 

-0, 

.961606 

16 

0. 

,304449 

1, 

,592619 

0. 

,146996 

0, 

,172964 

-0, 

.943442 

17 

0, 

,249485 

1, 

,650390 

0, 

,146996 

0, 

,214399 

-0, 

,920465 

18 

0, 

,198363 

1, 

,714437 

0. 

,146996 

0, 

,265094 

-0, 

.890942 

19 

0, 

,151840 

1, 

,786607 

0, 

,146996 

0, 

,328746 

-0, 

.851951 

20 

0, 

,110606 

1, 

,868882 

0. 

,146996 

0, 

,410545 

-0, 

.798852 

21 

0, 

,075269 

1, 

,965495 

0, 

,146996 

0, 

,519597 

-0, 

.722328 

22 

0, 

,046353 

2 

.082446 

0, 

,146996 

0, 

,668455 

-0, 

.603478 

23 

0, 

,024285 

2. 

,230950 

0, 

,146996 

0, 

,866926 

-0, 

.394782 

24 

0, 

,009388 

2, 

.419835 

0. 

,146996 

1, 

,009466 

0, 

.050859 

25 

0, 

,001884 

2. 

,339782 

0, 

,146996 

-0, 

,471799 

1, 

.214887 

26 

0, 

,001884 

-0. 

,156346 

0. 

,146996 

-4, 

.389847 

2, 

.320095 

27 

0, 

,009387 

-1, 

,322127 

0, 

,146996 

-3, 

,033496 

2, 

,003043 

28 

0, 

,024284 

-1, 

,560176 

0, 

,146996 

-2, 

,015200 

1. 

,727204 

29 

0. 

,046353 

-1, 

,622657 

0, 

,146996 

-1. 

,505832 

1. 

,569752 

30 

0, 

,075269 

-1 

,634560 

0. 

,146996 

-1 , 

,212790 

1 

.4-0549 

:: 

3. 

.  110606 

102 

0. 

,146996 

-1 , 

.021304 

. 

.401354 

.: 

0. 

..31340 

. 

.  6  x  Z  6  2  5 

3. 

,146996 

-0, 

,33E  • 

L 

.250007 

23 

0, 

.193363 

-1 , 

,596423 

0, 

,146996 

-0, 

,  730655 

— 

,309004 

34 

0, 

,249485 

-1. 

,578180 

0, 

,146996 

-0, 

,695048 

1, 

,274895 

35 

0, 

,304448 

-1. 

,560042 

0. 

,146996 

-0, 

,621833 

1. 

,245407 

36 

0, 

,362436 

-1, 

,542861 

0, 

,146996 

-0, 

,556431 

1. 

,218916 

37 

0, 

,422588 

-1. 

,526391 

0, 

,146996 

-0, 

,496616 

1 

,194569 

38 

0, 

,484005 

-1, 

.510876 

0, 

,146996 

-0, 

,440592 

1. 

.171595 

39 

0, 

,545766 

-1. 

.496317 

0, 

,146996 

-0, 

,387199 

1, 

.149437 

40 

0, 

,606941 

-1, 

.482638 

0, 

,146996 

-0, 

,335611 

1. 

.127617 

41 

0, 

,666606 

-1, 

.469497 

0, 

,146996 

-0, 

,285518 

1 

,105863 

42 

0. 

,723848 

-1 

.456874 

0, 

,146996 

-0, 

,236273 

1, 

.083745 

43 

0 

.777782 

-1 

.444335 

0, 

,146996 

-0, 

,187542 

1 

.060991 
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44 

45 
46 
47 
48 
49 
50 

0, 
0, 
0, 
0, 
0. 
0, 
0, 

,327561 
.372331 
.911495 
.944210 
.969875 
.987862 
,997429 

-1 
-1 

-1 
-X 
-1 
-1 

-1. 

.431955 
.419472 
.406547 
.393024 
.379863 
.363515 
.365097 

0.146996 
0.146996 
0.146996 
0.146996 
0.146996 
0.146996 
0.146996 

-0. 

-0 

-o. 

0 
0. 
0, 
0 

.138453 
.033061 
.035055 
.023027 
.091342 
.178655 
.303827 

1.037087 

1.011468 
0.933411 
0.951546 
0.912885 
0.361777 
0.734338 

CD 

- 

0.030956 

CL  = 

0.713821 

CM  =  - 

■0.190635 

TRAILING  VORTICES  DATA 


M 


X(M) 


Y(M) 


CIRC 


1 

2 

.383780 

0 

.232008 

-0 

.000933 

2 

2 

.333346 

0 

.225080 

-0 

.001625 

3 

2 

.234404 

0 

.216393 

-0 

.002229 

4 

2 

.235273 

0 

.206695 

-0 

.002784 

5 

2 

.136347 

0 

.i96301 

-0 

.003304 

6 

2 

.137600 

0 

.135376 

-0 

.003797 

7 

2 

.089004 

0 

.174067 

-0 

.004256 

8 

2 

.040442 

0 

.162550 

-0 

.004687 

9 

X 

.991957 

0 

.150853 

-0 

.005100 

10 

1, 

.943572 

0. 

.138989 

-0, 

.005481 

i  x 

1 , 

.395155 

0. 

.127193 

-0. 

.005301 

12 

t 

.  346653 

0, 

.115534 

-0, 

,006100 

13 

i! 

.798126 

0, 

,104170 

-0, 

,006349 

14 

i 

,749496 

0. 

.093078 

-0, 

,006559 

15 

i, 

.700793 

0, 

.082359 

-0. 

,006720 

16 

i, 

,651983 

0, 

.072087 

-0, 

,006339 

17 

1 

.603077 

0. 

.062343 

-0. 

,006896 

13 

x 

.554079 

0, 

.053196 

-0 

,006916 

19 

i! 

.505019 

0, 

.044648 

-0, 

.006385 

20 

l, 

.455917 

0, 

.036753 

-0. 

.006795 

21 

i, 

,406791 

0. 

,029591 

-0, 

,006643 

22 

i. 

,357635 

0, 

,023172 

-0 

,006443 

23 

l, 

,308651 

0. 

,017529 

-0, 

,006177 

24 

l. 

,259734 

0. 

,012684 

-0, 

,005852 

25 

i, 

,211021 

0, 

,008642 

-0, 

,005465 

26 

i , 

,162620 

0, 

,005401 

-0, 

,005014 

27 

l. 

,114706 

0, 

,002948 

-0, 

,004498 

28 

l, 

,067624 

0, 

,001210 

-0, 

,003917 

29 

l, 

,022530 

0, 

,000276 

-0, 

.003269 

TIME  STEP  TK  = 

1.999990 

TK  -  TKM1  ■ 

0.200000 

ALPHA(T)  =    7. 
U(T)      =    0. 

500000 
000000 

OMEGA (T) 
V(T) 

0.000000 
0.000000 

NITR 

VXW 

VYW 

WAKE 

THETA 

GAMMA 

0 

1 
2 

0.939783 
0.943367 
0.948647 

0.026143 

0.030223 
0.030031 

0.138029 

0.139770 
0.139825 

0.027811 

0.031363 
0.031699 

0.156187E+00 

0.160749E+00 
0.160765E+00 

CONVERGED  SOLUTION  OBTAINED  AFTER  NITR  =    2 

J 

X(J) 

Q(J) 

GAMMA 

CP(J) 

V(J) 

1 
2 

3 
4 
5 
6 
7 

0.997429 
0.987862 
0.969875 
0.944210 
0.911495 
0.872381 
0.827561 

1.116135 
1.110743 
1.126034 
1.150497 
1.179177 
1.210700 
1.244766 

0.160766 
0.160766 
0.160766 
0.160766 
0.160766 
0.160766 
0.160766 

0.320772 
0.221351 
0.158555 
0.116479 
0.086724 
0.065709 
0.051593 

-0.851401 
-0.907726 
-0.941336 
-0.962935 
-0.977640 
-0.987588 
-0.993866 
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3 

0 

.777784 

I 

.281123 

0 

.160766 

0 

.043212 

-0 

.997146 

9 

0 

.723850 

.319423 

0 

. 150766 

0 

.039694 

-0 

.997912 

10 

0. 

. 566609 

]    , 

.359943 

0, 

,160766 

0, 

,040813 

-0, 

,996296 

11 

0, 

.506945 

.402723 

0, 

,160766 

0, 

,046325 

-0, 

,992422 

-  -> 

0. 

,545769 

-J 

.448176 

0, 

,160766 

0, 

,056472 

-0, 

.986157 

13 

0, 

.434008 

i 

.496434 

0, 

,160766 

0, 

,071440 

-0, 

.977368 

14 

0, 

.422591 

i 

.547997 

0, 

,160766 

0, 

,091673 

-0, 

,965767 

15 

0. 

.362439 

i 

.503142 

0, 

,160766 

0, 

.117358 

-0, 

.950894 

+    .- 

0, 

,304449 

1 

.562397 

0, 

,160766 

0, 

.151404 

-0, 

.931843 

17 

0, 

.249435 

.727721 

0, 

,160766 

0, 

.193674 

-0, 

.907603 

13 

0, 

.198353 

i 

.793350 

0. 

,160766 

0, 

.247145 

-0, 

,876338 

19 

0, 

.151840 

i 
x 

.378113 

0, 

,160766 

0, 

,315652 

-0, 

.834967 

20 

0, 

.110606 

X   i 

.967481 

0, 

,160766 

0, 

,404369 

-0, 

.778575 

21 

0, 

,075259 

2 

.071099 

0, 

,160766 

0, 

,522025 

-0, 

.697329 

22 

0, 

.046353 

2 

.194821 

0, 

,160766 

0, 

,679633 

-0, 

.571292 

23 

0, 

,024235 

2 

.349161 

0. 

,  150766 

0, 

,381013 

-0, 

.350465 

24 

0, 

.009383 

2 

.539111 

0. 

,160766 

0, 

,937480 

0 

.119060 

25 

0, 

.001334 

2 

.420488 

0, 

,160766 

-0, 

,772853 

1 

.331881 

25 

0, 

.001884 

-0 

.236994 

0, 

,160765 

-4, 

,940613 

2 

.437122 

27 

0 

.009387 

-1 

.441366 

0. 

,160766 

-3, 

.294956 

2 

.071290 

28 

0, 

.024284 

-1 

.578374 

0, 

,160766 

-2, 

.145366 

i, 

.771569 

29 

0, 

.046353 

-1, 

.735029 

0, 

,160766 

-1, 

.575494 

1, 

.601988 

30 

0, 

.075269 

-1 

.740161 

0, 

,150766 

-1, 

,248166 

1 

.495590 

31 

0 

.110606 

-1 

.725694 

0, 

.160766 

-1 

.035354 

1 

.421368 

32 

0. 

.151340 

-1 

.705142 

0. 

,160766 

-0, 

.384674 

1 

.367020 

33 

0 

.  1 9  8  2  S  3 

-  ] 

.530851 

0, 

.160766 

-0, 

.770207 

1 

.323621 

34 

0 

.249435 

-1 

.655530 

0, 

.160766 

-0, 

.678837 

1 

.287744 

35 

0 

.304443 

-1 

.530341 

0 

.160766 

-0, 

.602871 

1 

.256978 

36 

0 

.  362436 

-1 

.606176 

0 

.160756 

-0, 

.537036 

1 

.229567 

37 

0 

.422538 

-1 

.582302 

0, 

.160766 

-0, 

.478581 

1 

.204600 

38 

0. 

.434005 

- 1 

.560513 

0 

.160766 

-0 

.425234 

1 

.181273 

39 

0 

.545766 

-1 

.539373 

0 

.160766 

-0 

.375344 

1 

.158993 

40 

0, 

.506941 

-I 

.519334 

0 

.160766 

-0 

.327562 

1 

.137220 

41 

0 

.  566605 

i 
i. 

.500277 

0, 

.150766 

-0, 

.231143 

1 

.115668 

42 

0 

.723848 

-1 

.432141 

0 

.160766 

-0, 

.235015 

1 

.093871 

43 

0, 

.777732 

-1 

.464664 

0, 

.160756 

-0 

.188536 

1 

.071527 

44 

0, 

.827561 

-1 

.448072 

0, 

.160766 

-0, 

.140512 

1 

.048044 

45 

0, 

.372381 

-1 

.432243 

0, 

.160766 

-0. 

.089803 

1 

.022803 

46 

0, 

.911495 

-1 

.417014 

0 

.160766 

-0, 

.035089 

0 

.995022 

47 

0, 

.944210 

-1 

.402392 

0, 

.  160766 

0 

.026050 

0 

.963233 

43 

0 

.969875 

.]_ 

.389331 

0, 

.160766 

0, 

.098463 

0 

.924430 

49 

0, 

.987362 

-1 

.379196 

0, 

,160766 

0 

.190332 

0 

.873000 

50 

0 

.997429 

-1 

.379105 

0 

.160766 

0, 

.319847 

0 

.795184 

CD 

~ 

0.022037 

CL  = 

1 

D. 709635 

CM  =  ■ 

■0. 

185327 

TRAILING  VORTICES  DATA 


M 

1 
2 
3 
4 
5 
6 
7 
3 

10 
11 
12 
13 
14 
15 
16 
17 
13 
19 
20 


2 

2 
2 
2 
2 
2 

2 

2 
2 

2 
2 
2 

2 
2 
2 
2 
2 


X(M) 

.923984 
.873369 
.823534 
.774192 
.725182 
.576492 
.523036 
.i79755 
.531597 
.433743 
.435886 
.387918 
.339993 
.291922 
.243772 
.195477 
.146968 
.098253 
.049394 
.000362 


0 
0 
0 
0 
0 
0 

0, 
0. 
0. 
0. 
0. 
0, 
0, 
0, 
0, 
0, 
0, 


Y(M) 

,318366 

,311733 

,302817 

,292439 

,281064 

,263830 

,256054 

,242831 

,229197 

215114 

200923 

186791 

172653 

158710 

144992 

131596 

113650 

106217 

094288 

082951 


CIRC 

0.000933 
0.001625 
0.002229 
0.002734 
0.003304 
003797 
004256 
004637 
.  05100 
0.005481 
0.005801 
0.006100 
0.006349 
0.006559 
0.006720 
0.006839 
0.006396 
0.006916 
0.006885 
0.006795 


■0 

■0 
•0 


1 


21 

1 

.951119 

0. 

.072314 

-0 

.006648 

22 

1 

.901708 

0 

.062412 

-0, 

.006443 

23 

1 

.852145 

0, 

.053298 

-0 

.006177 

24 

1 

.802427 

0 

.045045 

-0 

.005852 

25 

1 

.752599 

0 

.037677 

-0. 

.005465 

26 

1 

.702674 

0, 

.031252 

-0, 

.005014 

27 

1 

.652678 

0, 

.025824 

-0, 

.004498 

28 

1 

.602603 

0, 

.021481 

-0, 

.003917 

29 

1 

.552431 

0, 

.018390 

-0. 

.003269 

30 

1 

.501920 

0 

,017193 

-0, 

.002553 

31 

1 

.450588 

0, 

.016109 

-0, 

.002761 

32 

1. 

.378677 

0, 

.012411 

-0- 

.005504 

33 

1, 

.258439 

0, 

.007189 

-0, 

.007734 

34 

1 

.094864 

0, 

.003008 

-0, 

.009244 
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c.l       The  development  of  a 
computer  code  (U2DIIF) 
for  the  numerical  solu- 
tion of  unsteady,  invis- 
cid  and  imcompressible 
flow  over  an  airfoil. 


